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PREFACE 


The  work  described  in  this  report  was  authorized  under  Project  No.  0602622 A552. 
The  work  was  started  in  October  2000  and  completed  in  September  2001 . 

The  use  of  either  trade  or  manufacturers’  names  in  this  report  does  not  constitute  an 
official  endorsement  of  any  commercial  products.  This  report  may  not  be  cited  for  purposes  of 
advertisement. 


This  report  has  been  approved  for  public  release.  Registered  users  should  request 
additional  copies  from  die  Defense  Technical  Information  Center;  unregistered  users  should  direct  such 
requests  to  the  National  Technical  Information  Center. 


3 


Blank 


4 


CONTENTS 


1.  INTRODUCTION . 7 

2.  SMOKE  MECHANISMS  FOR  DEFEAT  OF  SENSORS . 7 

2. 1  Effective  Smoke  Screening:  Active  vs  Passive  and  Infrared  vs  Visible 

Systems . 7 

2.2  Smoke  Reductions  in  Signal  to  Noise  Ratio . 7 

2.3  Extinction  Coefficients . . 8 

3 .  CALCULATIONS  OF  EXTINCTION  COEFFICIENTS . 9 

3 . 1  Maximizing  Extinction  Coefficients  of  Polydispersions  of  Spheres . 9 

3.2  Calculating  Low  Frequency  Cross  Sections  for  Conductive  Rods  and  Disks . 11 

3.3  Calculating  High  and  Intermediate  Frequency  Extinction  Cross  Sections  for 

Conductive  Rods . 12 

3.4  Calculating  High  and  Intermediate  Frequency  Extinction  Cross  Sections  for 

Conductive  Disks . 15 

3 .5  Calculating  Conductivity  and  Permittivity  of  Metals  as  a  Function  of 

Frequency  Using  the  Drude  Model . 16 

3.6  Calculating  Size  Effects  on  Conductivity  of  Thin  Metal  Wires  and  Films . 16 

4.  DISCUSSION  OF  RESULTS  AND  FIGURES  FOR  DISKS  AND  RODS . 17 

4.1  Simple  Result  for  Some  Disks . 17 

4.2  Discussion  of  Size  Effects  on  Complex  Refractive  Index . 17 

4.3  Discussion  of  Rod  Infrared  Spectral  Extinction  Coefficients . 18 

4.4  Discussion  of  Disk  Infrared  Spectral  Extinction  Coefficients . 19 

4.5  Metal  Coated  Rods  and  Disks . 19 

4.6  Magnitude  of  Size  Effects  on  Extinction  Coefficients . 20 

4.7  Transition  Frequency  Errors  Near  the  Merger  of  the  High  and  Low 

Frequency  Solutions . 20 

5.  CONCLUSIONS . 20 

LITERATURE  CITED . 45 

SOURCE  CODE . 47 


5 


Blank 


6 


MAXIMIZING  INFRARED  EXTINCTION  COEFFICIENTS 
FOR  METAL  DISCS,  RODS,  AND  SPHERES 


1.  INTRODUCTION 

Infrared  spectral  extinction  coefficients,  extinction  cross  sections  per  volume  of 
material,  are  computed  for  metal  disks,  rods  and  spheres  as  a  function  of  dimensions, 
conductivity  and  wavelength.  Computations  for  rods  use  a  combination  of  infinite  cylinder  and 
Rayleigh  finite  and  perfectly  conducting  prolate  spheroid  scattering  theories.  Disk  computations 
use  a  combination  of  diffraction,  thin  film  optics,  Rayleigh  finite  and  perfectly  conducting  oblate 
spheroid  scattering  theories.  Sphere  computations  use  a  combination  of  Geometric  Optics, 
Anomalous  Diffraction,  and  Rayleigh  finite  and  perfect  conducting  sphere  scattering  theories. 
Complex  conductivities  are  predicted  as  a  function  of  wavelength  using  the  Drude  model  along 
with  thin  wire  and  thin  film  models  for  boundary  limited  electron  mean  free  path.  Spectral 
extinction  coefficients  are  maximized  using  extinction  coefficient  surface  plots  in  parameter 
hyperspace  to  identify  optimal  ranges  for  particle  dimensions  and  conductivities  at  a  particular 
wavelength.  Dramatic  differences  are  evident  in  the  way  that  essentially  one,  two  and  three 
dimensional  metal  particles  interact  with  infrared  radiation. 


2.  SMOKE  MECHANISMS  FOR  DEFEAT  OF  SENSORS 

2.1  Effective  Smoke  Screening;  Active  vs  Passive  and  Infrared  vs  Visible  Systems 

There  are  differences  between  smoke  screening  against  an  active  system 
compared  to  smoke  screening  against  a  passive  system  in  the  infrared.  Smoke  reduces  a  passive 
system  signal  (target  radiance)  by  one-way  attenuation  and  reduces  an  active  system  signal  by 
two-way  attenuation.  Smoke  clouds  superimpose  scattered  ambient  radiance  and  cloud  emitted 
radiance  onto  both  the  attenuated  signal  and  surrounding  attenuated  background  for  active  and 
passive  systems.  This  changes  the  noise  level,  which  is  proportional  to  the  square  root  of  signal 
and  background  radiance  for  ideal  photon  noise  limited  detection.  For  active  systems  the 
scattered  probe  radiance  is  removed  by  range  gating  and  scattered  ambient  radiance  and  cloud 
emitted  radiance  become  a  smaller  fraction  of  the  attenuated  signal  and  background  as  probe 
radiance  increases.  Unlike  the  visible  spectral  region  where  ambient  radiance  has  a  diffuse  sky 
component  representing  as  little  as  10%  of  the  incident  direct  solar  radiance,  ambient  infrared 
radiance  is  generally  more  isotropic  especially  under  rainy,  humid  or  cloudy  conditions  and 
therefore  less  influenced  by  differential  scatter  asymmetry  within  the  smoke.  Also  unlike  the 
visible  spectral  region  where  increasing  smoke  single  scatter  albedo  always  increases  the  smoke 
radiance  superimposed  onto  the  attenuated  signal  and  background,  in  the  infrared  the  combined 
effects  of  greater  scatter  and  less  absorption  will  result  in  less  cloud  radiance  when  the  drop  in 
thermal  emission  outweighs  the  increase  in  scattered  ambient  radiance.  This  occurs  when  smoke 
cloud  temperature  exceeds  ambient  temperatures. 

2.2  Smoke  Reductions  in  Signal  to  Noise  Ratio 

The  goal  of  effective  smoke  screening  is  to  reduce  signal  to  noise  ratio  using  the 
least  amount  of  smoke  material.  In  the  infrared  region  trade-offs  associated  with  increasing  noise 
by  adjusting  combinations  of  cloud  emittance,  reflectance  and  scattered  light  transmittance  versus 
reducing  signal  by  reducing  cloud  beam  transmittance  is  discussed  in  a  separate  report  (Embury, 
in  preparation).  In  the  visible  spectral  region  trade-offs  associated  with  increasing  noise  by 
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increasing  cloud  reflectance  and  scattered  light  transmittance  is  discussed  (Embury,  in 
preparation)  under  a  variety  of  solar  illumination  conditions  versus  reducing  signal  by  reducing 
cloud  beam  transmittance.  Those  reports  discuss  how  cloud  thermal  emission  increases  and  cloud 
reflectance  and  scattered  light  transmittance  decrease  as  smoke  particle  single  scatter  albedo 
decreases  and  the  signal  to  noise  ratio  tradeoffs  associated  with  changing  single  scatter  albedo, 
extinction  coefficient,  differential  scattering  pattern,  cloud  temperature  and  ambient  illumination. 

Minimizing  signal  by  reducing  cloud  beam  transmittance  using  a  minimal 
amount  of  smoke  material  is  discussed  here.  Signal  beam  transmittance  is  treated  by  breaking 
down  the  optical  depth  aCL  into  two  parts.  The  first  part  is  the  extinction  coefficient  a,  the 
extinction  (scatter  plus  absorption)  cross  section  per  volume  m2/cc  or  mass  m2/g  of  aerosol 
material.  The  second  part  is  the  concentration  path  length  product  CL,  the  volume  or  mass  of 
aerosol  per  cloud  area  in  the  image  plane,  where  C  is  the  aerosol  volume  cc/m3  or  mass  g/m3 
concentration  and  L  is  the  pathlength  in  meters  through  the  cloud.  Maximizing  the  extinction 
coefficient  minimizes  the  quantity  of  material  CL  required  to  reach  a  specified  optical  depth. 

2.3  Extinction  Coefficients 


Volume  extinction  coefficients  depend  on  optical  efficiency  factor  Q,  geometric 
cross  section  G  and  particle  volume  V  as  indicated  by  the  relationship 

a  =  QG/V. 

For  convex  particles  the  geometric  cross  section  averaged  over  random  orientation  is  one  fourth 
the  surface  area  S  so  that1 

«  =  QS/4V. 

Therefore,  the  volume  extinction  coefficient  can  be  increased  by  either  increasing  Q  or  by 
increasing  the  surface  area  to  volume  ratio.  In  order  to  significantly  increase  Q  over  a  narrow 
wavelength  range  we  can  choose  a  combination  of  particle  properties  that  produce  a  resonance. 
This  requires  particles  that  are  on  the  order  of  the  wavelength  in  size  and  that  fairly  monodisperse 
in  size  and  shape.  Increasing  the  surface  area  to  volume  ratio  means  reducing  particle  size  but 
eventually  the  particles  become  small  enough  that  increases  in  surface  to  volume  ratio  can  be 
offset  by  decreases  in  Q.  Spherical  particles  small  compared  to  the  wavelength  have  a  resonance 
at  m  =  (0,  V2)  which  broadens  to  a  wide  range  of  metallic  refractive  indices  for  rod  (fiber)  and 
disk  (flake)  shapes2.  This  resonance  maintains  Q  at  levels  around  2  in  the  case  of  metal  discs  and 
much  higher  levels  in  the  case  of  rods.  This  resonance  is  not  significantly  decreased  by  size 
distribution. 


Volume  extinction  coefficients  are  related  to  mass  extinction  coefficients  using 
particle  density  p 


a[m2/g]  =  a[m2/cc]//7[g/cc] 


and  it  is  convenient  to  express  the  surface  area  to  volume  ratio  in  micrometers  (i  to  obtain  the 
extinction  coefficient  directly  in  m2/cc.  It  should  be  mentioned  at  this  point  that  there  are  three 
ways  to  increase  the  mass  extinction  coefficient;  the  same  two  ways  that  were  just  mentioned  to 
increase  the  volume  extinction  coefficient  plus  reducing  the  particle  density  as  in  the  case  of 
hollow  or  foam  particles,  although  the  benefits  of  this  third  method  will  eventually  be  offset  by 
decreases  in  Q  at  very  low  densities. 
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The  remainder  of  this  report  involves  computation  of  the  volume  extinction 
coefficients  a  and  Q  which  are  functions  of  particle  size,  shape,  wavelength  and  complex 
refractive  index  which  in  turn  depends  on  constitutive  properties  as  well  as  wavelength  and 
minimum  dimension  in  the  case  of  nanoparticles. 


3.  CALCULATIONS  OF  EXTINCTION  COEFFICIENTS 

3.1  Maximizing  Extinction  Coefficients  of  Polvdispersions  of  Spheres 

Scattering  by  polydispersions  of  spheres  predicted  by  the  Mie  solution  can  be 
used  to  represent  scattering  by  collections  of  non-spherical  particles  as  long  as  they  are  roughly 
isometric  (equal  dimensions  along  any  three  orthogonal  axes)  and  have  size,  shape  and 
orientation  distributions  broad  enough  to  damp  out  resonances.  The  advantage  of  using  the 
sphere  solution  is  that  it  requires  less  computational  time  than  popular  non-spherical  techniques 
like  the  Extended  Boundary  Condition  Method,  the  Discrete  Dipole  Approximation  and  the 
Separation  of  Variables  Method  for  a  spheroid3. 

The  Mie  solution  for  the  scattering  cross  section  of  a  sphere  is1, 2 


2k_ 


Z(2«+i)(kl2 


and  the  Mie  solution  for  the  extinction  cross  section  of  a  sphere  is 


L(2»  +  l)Re(k,f 

■*;  K  n=l 


where  k  =  2k /X  and  for  a  sphere  with  complex  refractive  index  m  and  diameter  D  having  a  size 
parameter  x  =  KDjX , 


a  _  ”Wn (mx)y/ ( x )  -  y/n  (x)i/n  (mx) 

"  m¥n{mx)C(x)-Zn{x)¥'n{mx) 

b  _ 

Vn  (x)y/'n  (true) 

The  Riccati-Bessel  functions 

Vn(P)  =  Pin(P ').  tn(P)  =  Phn)(P) 


may  be  written  in  terms  of  the  Spherical  Bessel  functions 


Jn(P)  = 


(P) 
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which  are  defined  in  terms  of  the  Bessel  functions  of  the  first  kind  J  and  of  the  second  kind  Y. 


Spherical  aerosol  size  distributions  average  out  extinction  coefficient  resonances 
in  Q  yielding  a  simple  result  for  optimum  diameter  and  complex  refractive  index  producing 
maximum  achievable  extinction  coefficients  over  a  wavelength  band.  Mie  calculations  show  that 
the  maximum  achievable  extinction  coefficient  for  a  spherical  polydispersion  over  a  wavelength 
band  is  approximately 

«max[m2/cc]  =  3/Dopt[//], 

Dop,  is  the  minimum  diameter  maintaining  Q  >  2  over  the  wavelength  band.  Mie  calculations  also 
indicate  that  the  optimum  diameter  as  a  function  of  complex  refractive  index  m=(n,k)  and  the 
maximum  wavelength  in  the  band  is  approximately  the  same  result  obtained  using  the  theory 
of  Anomalous  Diffraction1  for  complex  refractive  indices  near  1  where  the  low  frequency  drop  in 
Q  occurs  around  a  phase  shift  of  n  through  the  sphere  diameter  yielding  the  following  result 

DoP1“^x/(2|ni-l|) 


This  indicates  that  in  the  limit  of  small  refractive  index  the  optimum  complex  refractive  index  is 
the  largest  achievable  since  that  minimizes  Dopt  and  maximizes  a.  For  refractive  indices  that  are 
not  small  Mie  calculations  show  that  Dopt  is  close  to  the  result  obtained  at  the  intersection  of  the 
low  frequency  Rayleigh  theory1 2  and  high  frequency  geometric  optics  plus  diffraction  scattering 
theories.1'2  Setting  the  high  frequency  extinction  cross  section  on  the  left  hand  side  of  the 
following  equation  equal  to  the  low  frequency  Rayleigh  absorption  plus  scatter  cross  sections  on 
the  right  had  side  we  get 


2k 


J 


f  2k  ^ 

3 

m2  —  1 1  8 

fz>f 

4k 

Im 

— = -  +-K 

— 

ItJ 

\  2  j 

_m2  +  2J  3  1 

U  J 

f kD Y 
\  A  y 


for  an  absorbing  and  scattering  sphere.  For  a  metal  with  little  absorption  we  drop  the  absorption 
term  (leading  term  on  the  right  hand  side  of  the  previous  equation)  and  add  a  magnetic  dipole 
scatter  term  equal  to  !4  the  electric  dipole  scatter  term  yielding  the  equation 


2k 


'  D ' 

2 

(S') 

8  (D'f 

-K  — 

3  V  2  J 

(kD  Y 

m2  - 1 

< 2 , 

V4  V 

l  X  J 

m2+  2 

Solving  this  equation  for  Dop,  we  have 


D 


Opt 


r3 

m2  +2 

2  > 

5 

V 

m2  —  1 

J 

A 

K 


Once  again  in  order  to  minimize  Dopt  to  maximize  a,  we  want  the  largest 
achievable  complex  refractive  index.  Dielectric  materials  such  as  titanium  dioxide  and  metals  at 
visible  and  longer  wavelengths  have  high  refractive  indices  and  achieve  maximum  spectral 


extinction  coefficients  of  ow  -  9/  A^.  Therefore  in  the  visible  spectral  region  where 
AmaxsO .7/i  we  have  c^x  =  13  m2/cc  and  in  the  infrared  where  Amax=14/r  we  have  a^x  **  0.64 
m2/cc  .  Conductive  foam  and  bubbles  can  increase  mass  extinction  coefficients  by  decreasing 
density  but  the  only  way  to  increase  both  mass  a[m2/g]  and  volume  a[m2/cc]  extinction 
coefficients  in  the  infrared  is  to  use  conductive  fibers  or  flakes. 

3.2  Calculating  Low  Frequency  Cross  Sections  for  Conductive  Rods  and  Disks 

Non-spherical  particle  scattering  theories4-20  such  as  the  Separation  of  Variables 
Method  for  spheroids,  the  Extended  Boundary  Condition  Method,  Discrete  Dipole 
Approximation  and  others  are  not  suited  to  compute  extinction  cross  sections  of  high  aspect  ratio 
prolate  and  oblate  spheroidal  shaped  particles  with  large  metallic  refractive  indices  because  of 
instabilities  in  the  solutions  and  long  computation  times.  Instead,  electromagnetic  extinction 
cross  sections  of  metal  disks  are  rapidly  computed  using  a  combination  of  diffraction,  thin  film 
optics  and  Rayleigh  oblate  finite  and  perfect  conducting  spheroid  scattering  theories. 
Electromagnetic  extinction  cross  sections  of  metal  rods  are  rapidly  computed  using  a  combination 
of  infinite  cylinder  and  Rayleigh  prolate  finite  and  perfect  conducting  spheroid  scattering 
theories.  These  theories  allow  computation  of  the  extinction  efficiency  factor  Q  for  polydisperse 
and  randomly  oriented  rods  and  disks  over  a  wide  wavelength  range.  Disk  and  rod  aerosol  size, 
aspect  ratio  and  orientation  distributions  average  out  extinction  cross  section  resonances 
permitting  disks  to  represent  lamellar  or  flake  shapes  and  permitting  rods  to  represent  filamentary 
or  fiber  shapes  leading  to  the  following  simplified  analysis.  As  mentioned  before,  Q  depends 
upon  particle  dimensions  and  complex  refractive  index  which  is  a  function  of  wavelength  and,  in 
the  case  of  very  small  particles,  also  a  function  of  particle  dimensions  which  limit  the  conduction 
electron  mean  free  path.  Complex  refractive  index  as  a  function  of  wavelength  is  predicted  using 
the  Drude  model  along  with  a  thin  film  model  for  boundary  reduction  of  conduction  electron 
mean  free  path  in  the  case  of  disks  and  along  with  a  wire  model  for  boundary  reduction  of 
conduction  electron  mean  free  path  in  the  case  of  rods.  In  the  analysis  the  orientation  distribution 
is  assumed  to  be  random  as  is  the  case  when  Brownian  and  turbulent  diffusion  forces  dominate 
over  preferred  orientation  forces  due  to  air  drag  acting  against  gravity. 

We  have  the  following  result  for  low  frequency  Rayleigh  absorption  and 
scattering  cross  sections  of  randomly  oriented  ellipsoids1'2  with  semi-axes  a,  b  and  c 

,  T  f|  1  1  1 

=  klm{-a,  +-<x  +- or, 

U  3  2  3  3 

(c”>=^lil“l|!+7K|!+7w!} 

2  71 

where  k  =  —  and  the  polarizability  cXj  along  semi-axis  i, 


a.  -  \7iabc- - — — — 


+  3Lj  (f ,  £m ) 


is  defined  in  terms  of  the  complex  permittivity  of  the  spheroid  £i,  that  of  the  medium  £m,  and  the 
depolarization  factor  L,  along  semi-axis  i. 
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The  complex  refractive  is  equal  to  the  square  root  of  the  complex  permittivity 

m  =  y/s  and  the  permittivity  of  air  is  close  to  unity  e„,  =  1.  For  a  prolate  spheroid  the 
depolarization  factor  along  the  symmetry  axis  is 


4 


l-e2 


r-l  +  —  In 
2e 


1  +6^ 
l-ej 


=  1- 


b2 


and  along  the  equatorial  axes  the  depolarization  factors  are 


4=(l-4)/2 


4  =  4 


For  an  oblate  spheroid  the  depolarization  factors  along  the  equatorial  axes  are 


4  = 


g2(*) 

2 


g(e)  = 


l-e 2 


4=4 


and  along  the  symmetry  axis  the  depolarization  factor  is 

4  =  1-24 

These  expressions  originate  from  the  Rayleigh  theory  of  scattering  by  ellipsoidal 
shapes  that  create  uniform  internal  electric  (polarization)  fields  in  a  uniform  external  electric 
field.  Prolate  spheroids  can  be  used  to  approximate  rods  and  the  oblate  spheroids  can 
approximate  disks.  This  theory  is  applicable  as  long  as  the  wavelength  both  inside  and  outside 
the  particle  is  large  compared  to  the  particle  dimensions.  Writing  rod  and  disk  major  dimension 
as  H  and  minor  dimension  as  h,  this  means  that  H  «  X  and  even  more  restrictively  that 
|mH| «  “K  in  the  case  of  metal  rods  and  disks  at  infrared  and  longer  wavelengths  where  m  can  be 
large.  We  often  expect  the  condition  on  mH  to  be  violated,  however  the  Rayleigh  electric  dipole 
results  above  stand  because,  unlike  spheres,  the  low  frequency  extinction  cross  section  correction 
due  to  the  magnetic  dipole  contribution  is  negligible  compared  to  the  electric  dipole  contribution 
for  rods  and  disks1. 

3.3  Calculating  High  and  Intermediate  Frequency  Extinction  Cross  Sections  for 

Conductive  Rods 


At  wavelengths  on  the  order  of  and  smaller  than  the  major  dimension  we  use  the 
infinite  cylinder  solution  for  rods  and  the  physical  optics  plus  diffraction  result  for  disks. 
Actually,  the  solution  giving  the  smallest  value  for  Q  is  chosen  at  wavelengths  above  and  below 
the  intersection  of  the  Rayleigh  theories  with  these  theories.  Either  Rayleigh  or  infinite  cylinder 
is  chosen  in  the  case  of  rods  and  either  Rayleigh  or  physical  optics  plus  diffraction  is  chosen  in 
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the  case  of  disks  on  either  side  of  the  intersection  of  these  two  theories  which  occurs  around 
X  «  3  H  If  a  switch  is  made  from  low  to  intermediate/high  frequency  solutions  at  X  »  3H,  there 
will  be  a  discontinuity.  This  is  fixed  with  a  bridging  formula 

(  V' 


a  = 


f  i  v 


V«L. 


V^IH  J 


where  n  □  1,  ot  is  the  low  frequency  and  am  is  the  intermediate/high  frequency  extinction 
coefficient.  I  have  chosen  simply  to  pick  the  smallest  of  the  two  at  any  wavelength  which  is 
equivalent  to  n  =  co  and  results  in  discontinuities  in  the  derivatives. 


The  solution  for  radiation  incident  upon  an  infinite  cylinder  at  an  oblique  angle 
of  incidence  is  known  but  I  have  chosen  to  use  the  computationally  less  time  consuming  result  for 
an  infinite  cylinder  illuminated  normal  to  its  axis  since  the  cross  section  for  an  orientation 
orthogonal  to  this  plane  is  negligible.  In  this  case  the  simple  closed  form  solution  for  the 
extinction  cross  section  for  unpolarized  radiation  incident  normal  to  the  cylinder  axis  can  be  used 
to  represent  the  extinction  cross  section  averaged  over  random  orientation  in  the  plane  normal  to 
the  Poynting  vector.  The  scatter  and  extinction  cross  sections  for  an  incident  electric  field 
polarized  in  the  plane  defined  by  the  incident  Poynting  vector  and  the  cylinder  axis  given  by2 


2sca,I 


^sca.  I  _  2 

2  aLI,  x 


u+2 


sK+w*) 

W— 1  f 


and  the  scatter  and  extinction  cross  sections  for  an  incident  electric  field  polarized  perpendicular 
to  that  plane  are 


Ssca  I] 


aon  f + 2S(Ku  r + \bM  r) 

n=l  7 


nil 


Qext,  II  ”  ]  ^011  +  2^0, 

X  L  n= 1 

and  for  unpolarized  incident  fields  we  have 


2sca  2  (^sctU  SscaJI  )»  Sex t  ^  (SexUI  Sext.1 ) 


where 
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For  computational  purposes  it  is  better  to  write  a„  and  b„  in  terms  of  the  log  derivative 

D..%e> 

JSp) 

noting  that  for  Bessel  functions  Z  =  Jor  Z=Y 


z:=z^(x)—z„(x) 

x 


thus 
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a  [fi,(^)/w  +  w/*R(*Wn-i(*) 

[Dn(mx)/m+n/x]H™(x)-H® ^(jc) 

6  [mDn(mx)+nlx]Jn{x)-Jn_x{x) 

[mDn{mx)  +  nlx\H^\x)-H^l{x) 

3.4  Calculating  High  and  Intermediate  Frequency  Extinction  Cross  Sections  for 

Conductive  Disks 


The  extinction  cross  section  of  a  thin  rectangular  plate  having  major  dimensions 
on  the  order  of  or  larger  than  the  wavelength  are  computed  based  on  Physical  Optics  and 
Diffraction  theories.  A  plate  of  area  A  illuminated  at  an  angle  9  with  respect  to  the  normal  to  the 
plate  surface  will  have  the  following  extinction  cross  section21 

Ce  =  2Acos0(\-T(0)) 

and  absorption  cross  section 

CA  =  Acos9(\-T{0)-R(9)) 

where  disk  reflectance  R{9)  and  transmittance  T(9)  are  those  for  an  infinite  slab22  illuminated  at 
an  angle  of  incidence  9 
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and  e  is  the  complex  permittivity  of  the  slab.  These  disk  cross  sections  must  be  integrated 
numerically  over  angle  to  obtain  the  cross  section  averaged  over  an  ensemble  of  randomly 
oriented  disks. 

3.5  Calculating  Conductivity  and  Permittivity  of  Metals  as  a  Function  of  Frequency 

Using  the  Prude  Model 


The  Drude  model  for  complex  refractive  index  can  be  used  with  reasonable 
accuracy  for  metals  at  infrared  and  longer  wavelengths23'27.  First  the  complex  conductivity  is 
written  as  function  of  DC  conductivity  cr,  frequency  oo,  and  conduction  electron  relaxation  time  t. 

a(co)  =  crl(\-i(OT) 

and  the  complex  permittivity  is  written 

e{co)  =  1  +  Ak  ia(co)  /  co 

so  that  the  complex  refractive  index  for  nonmagnetic  materials  becomes 
m(o)  =  e{co)V2 

3.6  Calculating  Size  Effects  on  Conductivity  of  Thin  Metal  Wires  and  Films 

The  effects  of  rod  boundaries  on  rod  complex  refractive  index  are  predicted  by  a 
theory28,29  for  thin  wire  conduction  electron  scatter.  Wire  boundaries  reduce  DC  conductivity  a 
from  its  bulk  value  c0  because  of  a  shorter  conduction  electron  relaxation  time  as  indicated  by  the 
following  expression 


^=i/0+i/*) 


where  K~~~  ls  rati°  °f  wire  diameter  to  conduction  electron  mean  free  path  A  in  the  bulk. 

The  effects  of  disk  boundaries  on  complex  refractive  index  is  given  by  a  theory28,30  for  thin  film 
conduction  electron  scatter.  Now  the  film  conductivity  ratio  becomes 


cr 


r5  K  K2  ^ 

v8  +  16-l6y 


where  /ris  the  ratio  of  film  thickness  to  conduction  electron  mean  free  path  in  the  bulk  and  £,  is 
the  exponential  integral 
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For  film  thicknesses  less  than  a  few  nanometers  quantum  corrections  may  be  necessary  but  these 
would  involve  experimental  measurements  to  determine  two  parameters  in  the  theory  for  each 
metal31,32. 


4.  DISCUSSION  OF  RESULTS  AND  FIGURES  FOR  DISKS  AND  RODS 

4.1  Simple  Result  for  Some  Disks 

For  thin  (h»H)  disks  opaque  to  the  radiation  with  diameters  greater  than  about 
one  third  the  wavelength  the  extinction  efficiency  factor  is  equal  to  two  and  the  surface  area  to 
volume  ratio  is  2/h  so  the  extinction  coefficient  may  be  written  simply  as 

«[m2/cc]  =  l/h  [  //] 

This  result  overestimates  the  extinction  cross  sections  of  disks  that  are  partially  transparent  and 
disks  with  major  dimensions  less  than  about  one  third  the  wavelength.  This  might  be  interpreted 
as  an  upper  limit  on  the  extinction  coefficient  of  a  disk  but  if  the  transmittance  T(9)  drops  off  less 
rapidly  than  l/h  it  is  not  actually  an  upper  limit. 

4.2  Discussion  of  Size  Effects  on  Complex  Refractive  Index 

Complex  refractive  indices  of  metals  in  bulk  form  where  particle  size  is  large 
enough  so  that  bounding  surfaces  have  negligible  effect  on  conduction  electron  mean  free  path 
and  relaxation  time,  have  been  computed  in  the  infrared  based  on  the  Drude  model  as  a  function 
of  wavelength  and  DC  conductivity.  Results  for  the  real  and  imaginary  parts  of  the  complex 
refractive  indices  appear  in  Figures  1  and  2  respectively  using  a  typical  value  of  r,/ao=2.5  ;  the 
radius  of  a  sphere  containing  one  conduction  electron  r,  in  units  of  the  Bohr  radius  ao ,  0.529  x 
10-8  cm.  This  then  sets  the  Fermi  velocity23 

Vf  =  4.2x1 08  /rsa0  [cm/sec] 

the  conduction  electron  relaxation  time 

r  =  0.22r  a03crx  1 (T14  [sec] 
and  the  conduction  electron  mean  free  path 

A  =  Vfr 

where  c  is  the  DC  conductivity.  When  specifying  DC  conductivity  a  [mho/cm]  we  must 
convert26  to  a  set  of  units  consistent  with  the  Drude  model  expressions  using 

0-[sec_1]  =  cr[mho/cm]c2[cm/sec]xl0‘9 

where  c  is  the  speed  of  light,  3  x  1010  cm/sec. 

Figures  1  and  2  show  the  real  and  imaginary  parts  of  metallic  refractive  indices  in 
the  bulk  form  as  a  function  of  wavelength  and  conductivity.  Rod  (thin  wire)  and  disk  (thin  film) 
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size  effects  on  complex  refractive  indices  are  shown  in  Figures  3-22  as  a  function  of  minimum 
dimension,  conductivity  and  wavelength.  Significant  differences  are  evident  between  the  infrared 
complex  refractive  indices  in  bulk  form  versus  both  thin  wire  and  film  forms  where,  even  up  to 
one  micrometer  minimum  dimension,  boundaries  reduce  conduction  electron  mean  free  path. 
Because  extinction  cross  sections  of  metal  rods  and  disks  depend  not  only  on  particle  major 
dimension,  minor  dimension,  and  wavelength  but  also  strongly  on  the  complex  refractive  index 
there  will  be  significant  differences  in  these  cross  sections. 

Particle  boundary  reductions  in  conduction  electron  mean  free  path  make  the 
complex  refractive  index  a  function  of  particle  size  as  well  as  conductivity  and  wavelength. 
Figures  3  and  4  show  significant  reductions  in  both  real  and  imaginary  parts  of  the  complex 
refractive  index  for  1  nm  diameter  wires  below  the  bulk  values  shown  in  the  previous  figures. 
Reductions  in  complex  refractive  indices  for  1  nm  thick  films  shown  in  Figures  5  and  6  are 
almost  as  great.  Figures  7  and  8  display  complex  refractive  indices  for  wires  with  a  conductivity 
of  6  x  105  mho/cm,  that  of  copper,  over  a  range  of  diameters  between  O.lnm  and  lOOOnm. 
Similarly  Figures  9  and  10,  1 1  and  12,  13  and  14  show  this  for  wires  at  conductivities  of  105 
mho/cm,  104  mho/cm  and  103  mho/cm  (that  of  graphite)  respectively.  At  the  higher  conductivities 
size  effects  are  apparent  up  to  wire  diameters  approaching  lOOnm.  Figures  15  and  16,  17  and  18, 
19  and  20,  21  and  22  show  complex  refractive  indices  for  films  with  conductivities  of  6  x  10s 
mho/cm,  105  mho/cm,  104  mho/cm  and  103  mho/cm  respectively  over  a  thickness  range  from 
O.lnm  to  lOOOnm.  Size  reductions  in  film  complex  refractive  indices  below  the  bulk  value  are 
evident  up  to  thicknesses  approaching  lOOnm  at  the  higher  conductivities. 

4.3  Discussion  of  Rod  Infrared  Spectral  Extinction  Coefficients 

Incorporating  the  size  effects  just  mentioned  into  calculations  of  spectral 
extinction  cross  sections  of  infinite  cylinders  and  Rayleigh  cylinders  make  it  possible  to 
accurately  predict  extinction  cross  sections  for  cylinders  with  diameters  less  than  lOOnm.  By 
combining  the  infinite  cylinder  solution  with  the  Rayleigh  solution  for  high  aspect  ratio  prolate 
spheroids  we  are  in  effect  combining  a  high  frequency  solution  for  rod  lengths  much  greater  than 
the  wavelength  with  a  low  frequency  solution  for  rod  lengths  much  less  than  the  wavelength.  In 
this  way  the  entire  extinction  cross  section  spectra  of  metal  rods  are  covered  as  long  as  extinction 
cross  section  resonances  in  the  region  where  rod  lengths  are  on  the  order  of  the  wavelength  get 
washed  out  by  averaging  length,  orientation  and  diameter  distributions.  This  should  occur  for 
sufficiently  broad  distributions  and  large  populations  of  metal  rods  as  would  be  found  for 
example  in  smoke  clouds.  Rod  lengths  H  and  diameters  h  appearing  in  the  figures  should  be 
interpreted  as  representative  of  the  length  and  diameter  population  distributions. 

Rod  extinction  coefficients  are  shown  in  Figures  23-68  as  a  function  of 
wavelength,  diameter,  length  and  conductivity.  Rod  length  and  conductivity  are  held  fixed  at 
levels  of  interest  in  Figures  23-38.  Rod  length  and  diameter  are  held  fixed  in  Figures  39-48, 
diameter  and  conductivity  are  held  fixed  in  Figures  49-60  and  wavelength  and  length  are  held 
fixed  in  Figures  61-68.  Extinction  coefficient  is  maximized  (~800m2/cc  at  a  wavelength  of  14 
micrometers)  at  an  optimum  combination  of  complex  refractive  index  that  is  the  largest 
achievable  (that  of  copper;  6xl05mho/cm),  any  length  greater  than  about  one  third  the  wavelength 
and  a  diameter  of  about  16nm.  Extinction  coefficients  nearly  that  large  (~600m2/cc  at  a 
wavelength  of  14  micrometers)  are  possible  for  metals  such  as  iron  with  a  conductivity  of 
105mho/cm  at  lengths  greater  than  about  one  third  the  wavelength  and  slightly  larger  diameters.  If 
lengths  that  large  are  not  achievable  then  extinction  coefficients  can  be  maximized  at  lower  levels 
using  optimum  combinations  of  length,  conductivity  and  diameter  at  various  wavelengths  read 
from  the  appropriate  figure. 
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4.4 


Discussion  of  Disk  Infrared  Spectral  Extinction  Coefficients 


We  assume  that  a  large  population  of  metal  disks  are  attenuating  infrared 
radiation  and  that  the  extinction  cross  section  is  averaged  over  sufficiently  broad  diameter, 
orientation  and  thickness  distributions  as  would  be  encountered  for  example  in  clouds  of  metal 
flake  smoke.  In  this  way  extinction  cross  section  resonances  are  wiped  out  in  the  region  where 
disk  diameter  is  on  the  order  of  the  wavelength  and  physical  optics  plus  diffraction  used  to 
compute  extinction  cross  sections  of  disks  with  diameters  large  compared  to  the  wavelength  is 
extended  up  to  wavelengths  on  the  order  of  the  diameter.  Likewise  the  Rayleigh  solution  for 
oblate  spheroid  scattering  used  to  compute  extinction  cross  sections  of  disks  with  diameters  small 
compared  to  the  wavelength  is  extended  down  to  wavelengths  on  the  order  of  the  diameter.  After 
incorporating  the  Drude  model  along  with  a  model  for  conduction  electron  mean  free  path 
limitations  due  to  thin  film  boundaries  the  full  extinction  cross  section  spectra  of  metal  disks  in 
the  infrared  are  covered  by  a  combination  of  these  high  and  low  frequency  theories.  Generally 
the  low  and  high  frequency  theories  are  found  to  meet  at  A  *  3D  .  Disk  diameters  H  and 
thicknesses  h  appearing  in  the  figures  should  be  interpreted  as  representative  of  the  diameter  and 
thickness  population  distributions. 

Disk  extinction  coefficients  are  shown  in  Figures  69-107  as  a  function  of 
wavelength,  diameter,  thickness  and  conductivity.  Disk  diameter  and  thickness  are  held  fixed  at 
levels  of  interest  in  Figures  69-77.  Disk  diameter  and  conductivity  are  held  fixed  in  Figures  78- 
93,  thickness  and  conductivity  are  held  fixed  in  Figures  94-99  and  wavelength  and  diameter  are 
held  fixed  in  Figures  100-107.  Extinction  coefficient  is  maximized  (~350m2/cc  at  a  wavelength 
of  14  micrometers)  at  an  optimum  combination  of  complex  refractive  index  that  is  the  largest 
achievable,  any  diameter  greater  than  about  one  third  the  wavelength  and  a  thickness  of  about 
lnm.  Extinction  coefficients  nearly  that  large  (~250m2/cc  at  a  wavelength  of  14  micrometers)  are 
possible  for  metals  such  as  iron  with  a  conductivity  of  105mho/cm  at  diameters  greater  than  about 
one  third  the  wavelength  and  slightly  larger  thicknesses.  If  diameters  that  large  are  not 
achievable  then  extinction  coefficients  can  be  maximized  at  lower  levels  using  optimum 
combinations  of  diameter,  conductivity  and  thickness  at  various  wavelengths  read  from  the 
appropriate  figure. 

4.5  Metal  Coated  Rods  and  Disks 


Rough  estimates  of  extinction  coefficients  of  metal  coated  dielectric  rods  and 
disks  can  be  made  by  reducing  the  previously  computed  extinction  coefficients  by  the  volume 
fraction  of  the  dielectric  substrate  material.  In  the  following  plots  this  was  done  assuming  a  metal 
rod  or  disk  equal  in  major  dimension  to  that  of  the  substrate  material  with  a  volume  equal  to  the 
coating  volume.  Metal  coated  dielectric  rod  extinction  coefficients  computed  in  this  way  are 
shown  in  Figures  108  and  109  as  a  function  of  coated  rod  overall  diameter  and  coating 
conductivity  at  a  wavelength  of  14  micrometers,  length  of  8  micrometers  and  dielectric  substrate 
diameter  of  20nm.  Metal  coated  dielectric  disk  extinction  coefficients  are  shown  in  Figures  1 10 
and  1 1 1  as  a  function  of  coated  disk  overall  thickness  and  coating  conductivity  at  a  wavelength  of 
14  micrometers,  diameter  of  7  micrometers  and  dielectric  substrate  thickness  of  20nm.  These 
figures  are  relevant  because  some  available  dielectric  flake  material  such  as  mica  and  talc  cleave 
preferentially  along  a  plane  that  further  reduces  thickness  during  milling.  Milling  and  then 
coating  these  particles  may  have  advantages  over  reducing  the  thickness  of  a  metal  flake  by  using 
its  malleability  during  milling  or  trying  to  produce  thin  free  standing  metal  flakes. 
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4.6 


Magnitude  of  Size  Effects  on  Extinction  Coefficients 


Metal  rod  extinction  coefficients  were  calculated  as  a  function  of  diameter  and 
length  at  a  conductivity  of  105mho/cm  and  a  wavelength  of  3|i  in  Figure  1 12  including  size 
effects  and  in  Figure  113  excluding  size  effects.  Large  decreases  in  the  extinction  coefficient  due 
to  size  effects  are  evident  for  diameters  under  a  few  nanometers.  A  similar  comparison  was  made 
at  a  wavelength  of  14p  in  Figures  1 14  and  115  and  again  extinction  coefficients  are  dramatically 
reduced  for  diameters  under  a  few  micrometers. 

Metal  disk  extinction  coefficients  were  calculated  as  a  function  of  thickness  and 
diameter  at  a  conductivity  of  105mho/cm  and  a  wavelength  of  3p  in  Figure  1 16  including  size 
effects  and  in  Figure  117  excluding  size  effects.  Large  decreases  in  the  extinction  coefficient  due 
to  size  effects  are  evident  for  thicknesses  under  a  few  nanometers.  A  similar  comparison  was 
made  at  a  wavelength  of  14p  in  Figures  1 18  and  119  and  again  extinction  coefficients  are 
dramatically  reduced  for  thicknesses  under  a  few  micrometers. 

4.7  Transition  Frequency  Errors  Near  the  Merger  of  the  High  and  Low  Frequency 

Solutions 


Low  and  high  frequency  extinction  cross  section  solutions  become  more  accurate 
the  farther  they  are  from  transition  or  intermediate  frequencies  located  where  the  solutions  are 
merged  using  one  of  several  bridging  techniques.  The  bridging  technique  used  throughout  this 
report  involves  choosing  the  smaller  extinction  coefficient  and  results  as  a  function  of  minimum 
and  maximum  dimensions  appear  in  Figures  1 12  and  1 14  for  rods  and  Figures  1 16  and  1 18  for 
disks.  This  approach  preserves  some  of  the  artificiality  of  the  construct  by  making  all  derivatives 
discontinuous  while  the  function  is  continuous. 

Comparable  results  involving  choosing  the  low  frequency  solution  for  major 
dimension  less  than  three  times  the  wavelength  appear  in  Figures  120  and  121  for  rods  and 
Figures  124  and  125  for  disks.  Here  the  function  itself  suffers  a  considerable  discontinuity  but 
the  high  and  low  frequency  extremes  are  not  contaminated  by  a  bridging  equation.  Finally  the 
reciprocal  bridging  equation  was  used  in  Figures  122  and  123  for  rods  and  Figures  126  and  127 
for  disks.  All  derivatives  are  continuous  and  any  evidence  of  the  bridge  has  vanished. 


5.  CONCLUSIONS 

Maximum  sphere  extinction  coefficients  in  units  of  m2/cm3  are  ~  9/X  where  X  is 
the  wavelength  of  interest  expressed  in  micrometers  (e.g.,  0.64  m2/cm3  at  X  =  14  micrometers). 
Optimum  diameter  and  complex  refractive  index  producing  this  maximum  is  a  diameter  D  ~  X/3 
and  the  largest  complex  refractive  index  achievable  typical  of  metals.  It  should  be  mentioned  that 

the  small  sphere  resonance  at  a  complex  refractive  index  of  gets  washed  out  over  any 

wavelength  band  just  as  the  large  real  refractive  index  resonances  get  washed  out  due  to  size 
distribution. 


The  most  dependable  metal  rod  and  disk  extinction  coefficient  calculations  are 
made  away  from  the  transition  region  where  major  dimension  is  on  the  order  of  the  wavelength. 
Accurate  low  frequency  results  appear  in  many  figures  where  the  major  dimension  is  less  than 
about  one  tenth  the  wavelength.  Accurate  high  frequency  results  appear  near  the  back  (major 
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dimension  much  greater  than  wavelength)  face  of  Figures  120,  121, 124  and  125  where  high 
frequency  extinction  coefficients  are  displayed  at  all  major  dimensions  greater  than  one  third  the 
wavelength  while  low  frequency  extinction  coefficients  are  displayed  at  major  dimensions  less 
than  one  third  the  wavelength.  These  four  figures  show  that  high  frequency  extinction 
coefficients  reach  maximum  values  of  ~300m2/cc  and  ~800m2/cc  at  wavelengths  of  3p  and  14p 
respectively  for  rods  and  ~150m2/cc  at  both  wavelengths  for  disks.  Optimum  rod  diameters  are 
~3nm  and  ~8nm  at  wavelengths  of  3p  and  14p  respectively  and  optimum  disk  thicknesses  are 
~0.7nm  at  both  wavelengths.  Higher  conductivity  metals  will  have  slightly  higher  maximum 
extinction  coefficients  and  smaller  optimum  minimum  dimensions  and  lower  conductivity  metals 
will  have  significantly  lower  maximum  extinction  coefficients  at  larger  optimum  minimum 
dimensions.  Rods  300nm  long  and  disks  300nm  in  diameter  have  maximum  low  frequency 
extinction  coefficients  that  are  almost  as  large  with  optimum  minimum  dimensions  close  to  those 
just  given  for  much  larger  major  dimension  rods  and  disks.  Further  reductions  in  rod  and  disk 
major  dimensions  reduce  the  maximum  extinction  coefficient  and  reduces  the  optimum  minimum 
dimension. 


There  are  essentially  three  types  of  errors  or  potential  errors  in  the  calculated 
extinction  cross  sections.  First  there  are  errors  that  result  from  extending  low  frequency 
calculations  outside  the  region  where  wavelength  is  much  greater  than  major  dimension  and 
extending  high  frequency  calculations  outside  the  region  where  wavelength  is  much  less  that 
major  dimension.  Unfortunately  a  great  deal  of  interest  is  directed  at  the  transition  region  where 
.major  dimension  is  on  the  order  of  the  wavelength  because  extinction  coefficients  increase  with 
increasing  major  dimension  until  major  dimension  is  on  the  order  of  the  wavelength.  Thus 
optimum  rod  length  and  optimum  disk  diameter  are  in  the  transition  region  and  no  bridging 
approach  can  completely  fix  this  problem.  Second  the  size  effects  models  are  accurate  down  to 
minimum  dimensions  around  lOnm  but  optimum  rod  diameters  are  in  the  3-8nm  range  and 
optimum  disk  thicknesses  are  less  than  lnm.  Third  the  calculations  assume  that  rods  have  a 
population  distribution  of  lengths  and  disks  have  a  distribution  of  diameters  that  average  out 
resonances  that  would  otherwise  occur  in  the  transition  region. 
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SOURCE  CODE 


Disk  and  Rod  Scattering  Codes 

Program  Darpadiskl 
0  =  310*10; 

lamb  =  lambmicron/1 0A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10*8/rsa0; 

A  =  H/h; 

taul  =  .22*mhopercml*10*-6*rsa0*3*10*-14; 
mfp  =  vftaul; 
kap  =  h/mfp; 
mhopercm  = 

mhopercml(l  -  .75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kap*2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul /mhopercm  1; 
sigO  =  10A-9  c*2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lamb)(l  -  Sin[th]A2/eps)A.5; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin(th]A2/eps)A.5  +  (1  -  Sin[th]*2/eps)*  5)/eps*  5; 
rr  =  I(zl  A2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)  10**4/01  +  hO),  (th,  10*-4,  Pi/2  -  10*-4}]/(Pi/2); 
ee  =  (1  -  1/A*2)*.5; 
ge  =  (l/ee*2  - 1)*.5; 

LI  =  ge/(2ee*2XPi/2  -  ArcTan[ge])  -  ge*2/2; 

L2  =  L1; 

L3  =  1  -  LI  -  L2; 

pi  =  (Pi/2)  H*2  (h  +  hOXeps  -  l)/(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  H*2  (h  +  hOXeps  -  l)/(3  +  3L2(eps  - 1)); 
p3  =  (Pi/2)  H*2  (h  +  hO)(eps  -  l)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

k*4(Abs[pl]*2  +  Abs[p2]*2  +  Abs[p3]*2)/(18  Pi))/((Pi/4)  H*2(h  + 
hO)  10*4); 

(♦lambmicron  =  10*ss;*) 
mhopercm  1  =  10*xx  ; 
h  =  3  10*-7; 

H  =  310*-5; 
h0  =  0; 
gl  =  Plot3D[ 

I f[ alphal  <=  alpha2,  alphal,  alpha2],  (lambmicron,  1, 14},  (xx,  3,  6}, 
PlotLabel  ->  "Metal  Disk  h=3nm  and  D=300nm", 

AxesLabel  ->  {"\[Lambda](\[Mu])M,  "log  \[Sigma]",  "\[  Alpha]  "}, 
PlotPoints  ->  20]; 
h  =  10*-6; 

H=  10*-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  (lambmicron,  1, 14},  (xx,  3, 6}, 
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PlotLabel  ->  "Metal  Disk  h=10nm  and  D=l\[Mu]M, 

AxesLabel  ->  {"\[Lambda]",  "log  \[Sigma]",  "\[  Alpha]  "},  PlotPoints  ->  20]; 
h  =  10^-6; 

H  =  210A-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  h=10nm  and  D=2\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "log  \[Sigma]",  "\[Alpha]  PlotPoints  ->  20]; 


Program  Darpadisk2 
c  =  3  10^0; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10A8/rsaO; 

A  =  H/h; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10A-14; 
mfp  =  vftaul ; 
kap  =  h/mfp; 
mhopercm  = 

mhopercml(l  -  ,75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul/mhopercml; 
sigO  =  10A-9  cA2  mhopercm; 
sig  =  sigO/  (1  -  I  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lambXl  -  Sinfth^/eps)^  5; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA.5; 
rr  =  I(zlA2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]>; 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =NIntegrate[ 

Cos(th](l  -  T)10^-4/(h  +  hO),  {th,  10M,  Pi/2  -  10A-4}J/(Pi/2); 
ee  =  (l-l/AA2)A.5; 
ge  =  (l/eeA2  - 1^.5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  L1; 

L3  =  1  -  LI  -  L2; 

pi  =  (Pi/2)  HA2  (h  +  h0)(eps  -  l)/(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  HA2  (h  +  h0)(eps  -  l)/(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  HA2  (h  +  h0)(eps  -  l)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kM(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(♦lambmicron  =  10Ass;*) 
mhopercml  =  10Axx ; 
h  =  1010A-7; 

H  =  5  10M; 
hO  =  0; 
gl  =  Plot3D[ 
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If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  h=10nm  and  D=5\[Mu]", 

AxesLabel  ->  {"\[Lambda](\[Mu])K,  'log  \[Sigma]",  "\[Alpha] 
PlotPoints  ->  20]; 

H  =  3  10^5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  h=10nm  and  D=300nm", 

AxesLabel  ->  {"\[Lambda](\[Mu])",  "log  \[Sigma]",  "\[Alpha] "}, 
PlotPoints  ->  20]; 

H  =  310**5; 

11  =  3  10**7; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  h=3nm  and  D=300nm", 

AxesLabel  ->  {"\[Lambda](\[Mu])",  "log  \[Sigma]",  ”\[Alpha] "}, 
PlotPoints  ->  20]; 

H  =  310A-5; 
h  =  10A-7; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3,  6}, 
PlotLabel  ->  "Metal  Disk  h=lnm  and  D=300nm", 

AxesLabel  ->  {"\[Lambda](\[Mu])",  "log  \[Sigma]",  "\tAlpha]  "}, 
PlotPoints  ->  20]; 

H=  10**5; 
h  =  3  10*-7; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  h=3nm  and  D=100nm", 

AxesLabel  ->  {"\[Lambda](\[Mu])",  "log  \[Sigma]",  "\[Alpha]  ”}, 
PlotPoints  ->  20]; 

H=  10^-5; 
h  =  10A-7; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3,  6}, 
PlotLabel ->  "Metal  Disk  h=lnm  and  D=100nm", 

AxesLabel  ->  {"\[Lambda](\[Mu])",  "log\[Sigma]",  ”\[Alpha]  "}, 
PlotPoints  ->  20]; 
h  =  3  10A-7; 

H=  10M; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel ->  "Metal  Disk  h=3nm  and  D=l\[Mu]", 

AxesLabel  ->  {"\[Lambda](\[Mu])M,  'log  \[Sigma]",  "\[Alpha]  "}, 
PlotPoints  ->  20]; 
h  =  3  10A-7; 

H  =  410a-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  h=3nm  and  D=4\[Mu]", 

AxesLabel  ->  {"V[Lambda](\[Mu])M,  "log  \[Sigma]'',  "\{Alpha]  ”}, 
PlotPoints  ->  20]; 
h  =  10  10A-7; 

H  =  410a-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
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PlotLabel  ->  "Metal  Diskh=10nm  and  D=4\[Mu]", 

AxesLabel  ->  {"\[Lambda](\[Mu])",  'log  \[Sigma]",  "\[Alpha] "}, 
PlotPoints  ->  20] 


Program  Darpadisk3 
c  =  3  10^10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=  4.210^8^530; 

A  =  H/h; 

taul  =  .22*mhopercm  1  *  1 0A-6*rsa0A3  *10^-14; 
mlp  =  vftaul; 
kap  =  h/mfp; 
mhopercm  = 

mhopercm  1(1  -  .75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul/mhopercm  1 ; 
sigO  =  10A-9  cA2  mhopercm; 
sig  =  sigO/  (1  -  I  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lambXl  -  Sintth^/eps)^; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA.5; 
rr  =  I(zlA2  -  z2A2)T  an[a2]/(2z  1  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)10'M/(h  +  hO),  (th,  10^,  Pi/2  -  10A-4}]/(Pi/2); 
ee  =  (l-l/AA2)A.5; 
ge  =  (l/eeA2  -  iy\5; 

LI  =  ge/(2eeA2XPi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  L1; 

L3  =  l  -LI  -L2; 

pi  =  (Pi/2)  HA2  (h  +  hO)(eps - 1)/(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  HA2  (h  +  hO)(eps  -  l)/(3  +  3L2(eps  - 1)); 
p3  =  (Pi/2)  HA2  (h  +  hO)(eps  - 1)/(3  +  3L3(eps  - 1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kM(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(♦lambmicron  =  lO^s;*) 
h  =  10Axx; 

mhopercm  1  -  6  10A5; 

H  =  4  10A-4; 
h0  =  0; 
gl  =  Plot3D[ 

Ifjalphal  <=  alpha2,  alphal,  alpha2],  (lambmicron,  1,  14},  (xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=6xlOA5  and  D=4\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[Alpha]  PlotPoints  ->  20]; 
mhopercm  1  =  10A5; 
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gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmiCTon,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10A5  and  D=4\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[Alpha] "},  PlotPoints  ->  20]; 
mhopercml  =  10*4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {laxnbmiCTon,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10A4  and  D=4\[Mu]", 

AxesLabel  ->  ("\[Lambda]",  "logh",  "\[Alpha] "},  PlotPoints  -> 20]; 
mhopercml  =  10*3; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmiCTon,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10*3  and  D=4\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[Alpha] "},  PlotPoints  ->  20]; 
mhopercml  =  6  10*5; 

H  =  10*-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=6xl0A5  and  D=l\[Mu]H, 

AxesLabel  ->  {"\[Lambda]n,  "logh",  "\[ Alpha] "},  PlotPoints  -> 20]; 
mhopercml  =  10*5; 

H=10*-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10*5  and  D=l\[Mu]", 

AxesLabel  ->  {^[Lambda]",  "logh",  ”\[Alpha]  "},  PlotPoints  -> 20]; 
mhopercml  =  10*4; 

H  =  10*-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10*4  and  D=l\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "logh",  ”\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10*3; 

H  =  10*-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10*3  and  D=l\[Mu]", 

AxesLabel  ->  {^[Lambda]",  "logh",  "\[Alpha] "},  PlotPoints  ->  20] 


Program  Darpadisk4 
c  =  3  10*10; 

lamb  =  lambmicron/1 0*4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10*8/rsa0; 

A  =  H/h; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10*-14; 
mfp  =  vftaul ; 
kap  =  h/mfp; 
mhopercm  = 

mhopercml(l  - .75(kap  -  kap*3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kap*2/16)Exp[-kap]) ; 
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tau  =  mhopercm  taul/mhopercm  1 ; 

sigO  =  10^-9  mhopercm; 

sig  =  sigO/  (1  - 1  omeg  tau); 

eps  =  1  + 1 4  Pi  sig/omeg; 

a 2  =  (2Pi  h  epsA.5/lamb)(l  -  Sin[th]A2/eps)A  5; 

zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5yepsA.5; 
rr  =  I(zlA2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alpha  1  =NIntegrate[ 

Cos[th](l  -  TJIOM/Oi  +  hO),  {th,  10M,  Pi/2  - 
ee  =  (1  -  l/A'S)''.  5; 
ge  =  (l/ee^2  -1^.5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  L1; 

L3  =  1  -  LI  -L2; 

pi  =  (Pi/2)  HA2  (h  +  hOXeps  - 1)/(3  +  3Ll(eps  -  1)); 

p2  =  (Pi/2)  HA2  (h  +  hOXeps  - 1)/(3  +  3L2(eps  -  1)); 

p3  =  (Pi/2)  HA2  (h  +  hO)(eps  -  l)/(3  +  3L3(eps  -  1)); 

alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

k^Absfpl]^  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(♦lambmicron  =  10Ass;*) 
h  =  10Axx; 

mhopercm  1  =  6  10^; 

H  =  3  10A-5; 
hO  =  0; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alpha  1,  alpha2],  {lambmicron,  1, 14},  (xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=6xlOA5  and  D=300nm", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[ Alpha]  PlotPoints  ->  20]; 

mhopercm  1  =  10^5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  ^Sigma^lO^  and  D=300nm", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercm  1  =  10A4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10A4  and  D=300nm", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[Alpha] "},  PlotPoints  ->  20]; 
mhopercm  1  =  10A3; 
gl  =  Plot3D[ 

Ifjalphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=l(Kv3  and  D=300nm", 

AxesLabel  ->  {"\[Lambda]",  "logh",  ”\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercm  1  =  6  10^; 

H=  1.00001  10^5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=6xlOA5  and  D=100nm", 

AxesLabel  ->  {"\[Lambda]",  "logh”,  ”\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercm  1  =  10^5; 
gl  =  Plot3D[ 
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Iffalphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10A5  and  D=1  OOnm", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[Alpha] PlotPoints  ->  20]; 
mhopercml  =  10*4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  (xx,  -5,  -8}, 
PlotLabel ->  "Metal  Disk  \[Sigma]=10A4  and  D=1  OOnm", 

AxesLabel  ->  {"\[Lambda]",  "logh",  '"{Alpha]  "},  PlotPoints  20]; 
mhopercml  =  10*3; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Disk  ^SigmaHO^  and  D=100nm", 

AxesLabel  ->  {"\[Lambda]",  "logh",  "\[Alpha]  "},  PlotPoints  -> 20] 


Program  Darpadisk5 
c  =  310A10; 

lamb  =  lambmicron/ 1 0A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10A8/rsa0; 

A  =  H/h; 

taul  =  .22*mhopercm  1  *  1 0A-6*rsaCK'3  *  1 0*- 1 4; 
mip  =  vPtaul; 
kap=h/mfp; 
mhopercm  = 

mhopercml  (1  -  .75(kap  -  kapA3/l 2)ExpIntegralEi [-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul /mhopercml; 
sigO  =  10A-9  c*2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lambXl  -  Sin[th]A2/eps)A.5; 
zl  =  (1/Costth]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA.5; 
rr  =  I(zlA2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R=  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)10M/(li  +  h0),  {th,  10M,  Pi/2  -  10A-4}]/(Pi/2); 
ee  =  (l-l/AA2)A.5; 
ge  =  (l/eeA2  - 1^.5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  LI; 

L3  =  1  -  LI  -  L2; 

pi  =  (Pi/2)  HA2  (h  +  hOXeps  -  l)/(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  HA2  (h  +  hOXeps  -  l)/(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  HA2  (h  +  hOXeps  - 1)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Imfpl  +  p2  +  p3]/3  + 

kA4(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(♦lambmicron  =  10Ass;*) 
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h  =  10**7; 

mhopercm  1  -  6  10*5; 

H=  1.00001  10*xx; 
h0  =  0; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=6xlO*5  and  h=lnm", 

AxesLabel  ->  {"\[Lambda]",  "log  D",  "\[Alpha]  "},  PlotPoints  ->  20] 


Program  Darpadisk6 
c  =  3  10*10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10*8/rsa0; 

A  =  H/h; 

taul  =  .22*mhopercml*10*-6*rsa0*3*10*-14; 
mfp  =  vftaul ; 
kap  =  h/mfp; 
mhopercm  = 

mhopercm  1(1  -  ,75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kap*2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul/mhopercml; 
sigO  =  10A-9  c*2  mhopercm; 
sig  =  sigO/  (1  -  I  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lamb)(l  -  Sin[th]*2/eps)*.5; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA.5; 
rr  =  I(zl A2  -  z2A2)T an[a2]/(2z  1  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zl A2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)10A-4/(h  +  hO),  {th,  10**4,  Pi/2  -  1  0*-4}]/(Pi/2); 
ee  =  (l  -  1/A*2)*.5; 
ge  =  (1/66*2  - 1)*.5; 

LI  =  ge/(2ee*2)(Pi/2  -  ArcTanfge])  -  ge*2/2; 

L2  =  L1; 

L3  =  l  -LI  -L2; 

pi  =  (Pi/2)  H*2  (h  +  hO)(eps  - 1)/(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  H*2  (h  +  hO)(eps  -  l)/(3  +  3L2(eps  - 1)); 
p3  =  (Pi/2)  H*2  (h  +  hOXeps  - 1)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

k*4(Abs[pl]*2  +  Abs[p2]*2  +  Abs[p3]A2)/(l  8  Pi))/((Pi/4)  H*2(h  + 
hO)  10*4); 

(♦lambmicron  =  10*ss;*) 
h  =  10  10*-7; 
mhopercm  1  =  6  10*5; 

H=  1.00001  10*xx; 
hO  =  0; 
gl  =  Plot3D[ 
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If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Disk  VtSigmaHxlO^  and  h=10nm", 

AxesLabel  ->  {"\[Lambda]",  "log  D",  "\[Alpha]  "},  PlotPoints  ->  20]; 
h  =  1(^-7; 
mhopercml  =  10A3; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  (lambmicron,  1, 14},  (xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=lCK'3  and  h=lnm", 

AxesLabel  ->  ("\[Lambda]",  "log  D",  "\[Alpha]  "},  PlotPoints  ->  20]; 
Show[%,  PlotRange  ->  (0,  50}]; 

Show[%,  PlotRange  ->  (0, 15}]; 

Show[%,  PlotRange  ->  (0,  20}]; 
h=  10  10*-7; 
mhopercml  =  10*3; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  (lambmicron,  1, 14},  (xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Disk  VfSigmablO^  and  h=10nm", 

AxesLabel  ->  {^[Lambda]",  "log  D",  "\[Alpha]  "},  PlotPoints  ->  20]; 
Show[%,  PlotRange  ->  (0,  20}]; 
h=  10  10A-7; 
mhopercml  =  10A5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  (lambmicron,  1, 14},  (xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10A5  and  h=10nm", 

AxesLabel  ->  ("\[Lambda]",  "log  D",  "\[ Alpha]  "},  PlotPoints  ->  20]; 
h  =  10 10A-7; 
mhopercml  =  10^4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  (lambmicron,  1, 14},  (xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10A4  and  h=10nm", 

AxesLabel  ->  {^[Lambda]",  "log  D”,  "\[Alpha]  "},  PlotPoints  ->  20]; 
Show[%,  PlotRange  ->  (0,  50}] 


Program  Darpadisk7 
c  =  3  10*10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.210A8/rsa0; 

A  =  H/h; 

taul  =  .22*mhopercm  1  *10A-6*rsa0A3*  1 0A-14; 
mfp  =  vPtaul ; 
kap  =  h/mfp; 
mhopercm  = 

mhopercml(l  -  .75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul /mhopercml; 
sigO  =  10^-9  cA2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lamb)(l  -  Si^th^/eps)^; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 
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z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.  5)/epsA.5; 
rr  =  I(zl  A2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alpha  1  =  NIntegrate[ 

Cos[th](l  -  T)1  (^-4/(11  +  hO),  {th,  10^-4,  Pi/2  -  10A-4}]/(Pi/2); 
ee  =  (1  -  l/A*2y-.5; 
ge  =  (1^2- 1^.5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTanfge])  -  geA2/2; 

L2  =  LI; 

L3  =  1  -  LI  -  L2; 

pi  =  (Pi/2)  HA2  (h  +  hO)(eps  - 1)/(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  HA2  (h  +  hO)(eps  - 1)/(3  +  3L2(eps  - 1)); 
p3  =  (Pi/2)  HA2  (h  +  hO)(eps  - 1)/(3  +  3L3(eps  - 1)); 
alpha2  =  (k  Imjpl  +  p2  +  p3]/3  + 

k^Abstpl]^  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(♦lambmicron  =  10Ass;*) 
mhopercml  =  10Axx ; 
lambmicron  =  14; 
h  =  10Ayy; 

H  =  4  1(^-4; 
h0  =  0; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=14\[Mu]  and  D=4\[Mu]", 
AxesLabel  ->  {"logh",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 
H=10A-4; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=14\[Mu]  and  D=l\[Mu]", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 
H  =  3  10A-5; 

gl  =  Plot3D[lf[alphal  <=  alpha2,  alphal ,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=14\[Mu]  and  D=300nm", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 
H=10A-5; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal ,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=14\[Mu]  and  D=100nm", 
AxesLabel  ->  {"log  h",  'log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20] 


Program  Darpadisk8 
c  =  3  10A10; 

lamb  =  lam  bm  i  cron/ 1 0A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf  =  4.2  1 0A8/rsaO; 

A  =  H/h; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10A-14; 
mfp  =  vftaul; 
kap  =  h/mfp; 
mhopercm  = 

mhopercml(l  -  .75(kap  -  kapA3/I2)Exp!ntegralEi[-kap]  - 
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3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul/mhopercm  1 ; 
sigO  =  10^-9  c*2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lamb)(l  -  Sinfth^/eps^.S; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA  5; 
rr  =  I(zl A2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zl A2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)10A-4/(h  +  hO),  {th,  10M,  Pi/2  -  10M}]/(Pi/2); 
ee  =  (l-l/AA2)A.5; 
ge  =  (l/eeA2  -  iy\5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  LI; 

L3  =  1  -  LI  -  L2; 

pl  =  (Pi/2)  HA2  (h  +  hOXeps  -  l)/(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  HA2  (h  +  hO)(eps  - 1)/(3  +  3L2(eps  - 1)); 
p3  =  (Pi/2)  HA2  (h  +  hO)(eps  - 1)/(3  +  3L3(eps  - 1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kA4(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(•lambmicron  =  KK'ss;*) 
mhopercm  1  =  KK'xx ; 
lambmicron  =  14; 
h=  lO'^yy; 

H=  7  lOM; 
h0  =  0; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5),  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=14\[Mu]  and  D=7\[Mu]", 
AxesLabel  ->  {"logh",  "log  \[Sigma]",  ”\[Alpha]"},  PlotPoints  ->  20]; 
gl  =  Plot3D[alphal,  {yy,  -8,  -5},  {xx,  3, 6}, 

PlotLabel  ->  "Metal  Disk  \[Lambda]=14\[Mu]  and  D=7\[Mu]", 
AxesLabel  ->  {"logh",  "log \[Sigma]”,  "\[ Alpha]"},  PlotPoints -> 20]; 
h0  =  20  10A-7; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Coated  Disk  \(Lambda]=14\[Mu]di0=20nmJ>=7\[Mu]", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20] 


Program  Darpadisk9 
c  =  3  10A10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf  =  4.2  lO^/rsaO; 

A  =  H/h; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10A-14; 
mfp  =  vPtaul; 
kap  =  h/mfp; 
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mhopercm  = 

mhopercm  1(1  -  .75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul/mhopercm  1 ; 
sigO  =  10^9  c^2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lamb)(l  -  Sin[th]A2/eps)A.5; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA.5; 
rr  =  I(zl A2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)1 0A-4/(h  +  hO),  {th,  1(^-4,  Pi/2  -  10^-4}]/(Pi/2); 
ee  =  (1  -  l/A^.  5; 
ge  =  (l/eeA2-l)A.5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  L1; 

L3  =  l  -LI  -L2; 

pi  =  (Pi/2)  HA2  (h  +  hO)(eps  - 1)/(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  HA2  (h  +  hO)(eps  -  l)/(3  +  3L2(eps  -  1 )); 
p3  =  (Pi/2)  HA2  (h  +  hO)(eps  -  l)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kM(Abs[p  1  ]A2  +  Abs[p2] A2  +  Abs[p3]A2)/(  1 8  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(•lambmicron  =  1  O^ss;*) 
mhopercm  1  =  10Axx  ; 
lambmicron  =  14; 
h  =  10Ayy; 

H  =  5  1(^-4; 
h0  =  0; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5),  (xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=14\[Mu]  and  D=5\[Mu]", 
AxesLabel  ->  ("log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 
lambmicron  =  7; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=7\[Mu]  and  D=5\[Mu]", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]H,  "\[Alpha]"},  PlotPoints  ->  20]; 
lambmicron  =  5; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=5\[Mu]  and  D=5\[Mu]", 
AxesLabel  ->  {"log  h”,  "log  \[Sigma]",  ”\[Alpha]"},  PlotPoints  ->  20]; 
lambmicron  =  3; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Disk  \[Lambda]=3\[Mu]  and  D=5\[Mu]”, 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 


Program  DarpadisklO 
c  =  3  10A10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
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k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=  4.2  lO^/rsaO; 

A  =  H/h; 

taul  =  .22*mhopCTcml*10A-6*rsa0A3*10A-14; 
mfp  =  vftaul ; 
kap  =  h/mip; 
mhopercm  = 

mhopercml(l  -  .75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap])  ; 
tau  =  mhopercm  taul/mhopercml; 
sigO  =  10^9  <^2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lambXl  -  Smph^/epsr.S; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA.5; 
rr  =  I(zl  A2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zlA2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)10A-4/(h  +  h0),  {th,  lO'M,  Pi/2  -  10A-4}]/(Pi/2); 
ee  =  (l-l/AA2)A.5; 
ge  =  (l/eeA2  - iy\ 5; 

LI  =  ge/(2eeA2XPi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  L1; 

L3  =  1  -  LI  -  L2; 

pi  =  (Pi/2)  HA2  (h  +  hOXeps  -  l)/(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  HA2  (h  +  hOXeps  -  l)/(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  HA2  (h  +  hOXeps  -  l)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kA4(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  HA2(h  + 
hO)  10A4); 

(*lambmicron  =  lO^s;*) 
mhopercm  1  =  10Axx ; 
lambmicron  =  14; 
h=10Ayy; 

H  =  7  1(^-4; 
h0=1010A-7; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Coated  Disk  \[Lambda]=14\[Mu]Ji0=10nm4>7\[Mu]", 
AxesLabel  ->  {"log  h",  'log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20] 


Program  Darparodla 
0=310^0; 

lamb  =  lambmicron/1 0A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10A8/rsa0; 

A  =  H/h; 
x  =  kh; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10A-14; 
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mfp  =  vPtaul ; 
kap  =  h/mfp; 

mhopercm  =  mhopercml/(l  +  mfp/h); 
tau  =  mhopercm  taul /mhopercm  1 ; 
sigO  =  10**9  c*2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
m  =  eps*0.5; 
z  =  m  x; 

qel  =  (2/x)Re[((m  Bessel  J[-l,  z]/BesselJtO,  z])  Bessel J[0,  x]  - 
BesselJ[-l, 

x])/((m  BesselJ[-l,  z]/BesselJ[0,  z]XBesselJ[0,  x]  + 

I  Bessel Y[0,  x])  -  (BesselJ[-l,  x]  + 

I  Bessel Y[-l,x]))  + 

2Sum[((m  Bessel J[j  - 1,  z]/BesselJ[j,  z])BesselJO,  x]  - 
Bessel  J[j  - 1, 

x])/((m  BesselJ[j  - 1,  z]/BesselJ[j,  z])(BesselJ[j,  x]  + 

I  Bessel  Y[j,  x])  -  (Bessel J[j  - 1,  x]  + 

I  BesselY[j  - 1,  x])),  (j,  1, 19}]]; 
qe2  =  (2/x)Re[(BesseIJ[-l,  z]  BesselJ[0,  x]/(m  BesselJ[0,  z])  - 
BesselJ[-l, 

x])/(  (BesselJ[-l,  z]/(m  BesselJ[0,  z]))(BesselJ[0,  x]  + 

I  BesselY[0,  x])  -  (BesselJ[-l ,  x]  + 

I  BesselY[-l,  x]))  + 

2Sum[((  BesselJ[j  -  1,  z]/(m  BesselJ[j,  z])  -  j  /(m  z)  + 
j/x)BesselJ[j,  x]  - 
Bessel J[j  - 1, 

x])/(  (Bessel J[j  - 1,  z]/(m  BesselJ[j,  z])  -  j/(m  z)  + 
j/x)(BesselJ[j,  x]  +  I  BesselY[j,  x])  -  (BesselJ[ 
j  - 1,  x]  + 1  BesselYjj  *  1,  x])),  (j,  1, 19}]]; 
qe  =  (qel  +  qe2)/3; 
alphal  =  qe/(Pi  h  10*4/4); 
ee  =  (1  -  1/A*2)*.5; 

LI  =  (l/ee*2  -  1X-1  +  Log[(l  +  ee)/(l  *  ee)]/(2ee)); 

L2  =  (1  -  LI  )/2; 

L3  =  L2; 

pi  =  (Pi/2)  H  h*2(eps  -  iy(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  H  h*2(eps  -  iy(3  +  3L2(eps  -  1)); 

P3  =  (Pi/2)  H  h*2(eps  -  iy(3  +  3L3(eps  *  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

k*4(Abs[pl]*2  +  Abs[p2]*2  +  Abs[p3]*2)/(18  Pi)y((Pi/ 

4)  Hh*2  10*4); 
h  =  30  10**7; 

H  =  4  10**4; 
mhopercm  1  =  10*xx ; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=4\[Mu]  and  h=30nm", 

AxesLabel  ->  {"\[Lambda]M,  "log  \[Sigma]",  ”\[Alpha]  "}, 

PlotPoints  ->  20]; 
h  =  10  10**7; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1,  14},  {xx,  3, 6}, 
PlotLabel ->  "Metal  Rod  H=4\[Mu]  and  h=10nm", 

AxesLabel  ->  {"\[Lambda]",  "log  \[Sigma]”,  "\[Alpha]  "}, 

PlotPoints  ->  20]; 
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h  =  10  10**7; 

H  =  10*-4; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=l\[Mu]  and  h=10nm", 

AxesLabel  ->  {"\[Lambda]",  "log  \[Sigma]",  "\[Alpha]  ■}, 

PlotPoints  ->  20]; 
h  =  3010A-7; 

H=  10A-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  (xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=l\[Mu]  and  h=30nm", 

AxesLabel  ->  {”\[Lambda]",  'log  \[Sigma]",  "\[Alpha]  "}, 

PlotPoints  ->  20]; 
h  =  30 10A-7; 

H  =  3  10A-5; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=300nm  and  h=30nm", 

AxesLabel  ->  {^[Lambda]",  'log  \tSigma]",  "\[Alpha]  "}, 

PlotPoints  ->  20]; 
h=1010A-7; 

H  =  310*-5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=300nm  andh=10nm", 

AxesLabel  ->  {^[Lambda]",  'log  \[Sigma]",  "\[Alpha]  "}, 

PlotPoints  ->  20]; 
h  =  10A-7; 

H  =  310A-5; 
gl  =  Plot3D[ 

I f[ alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=300nm  andh=lnm", 

AxesLabel  ->  {^[Lambda]",  'log  \[Sigma]",  "\{Alpha]  "}, 

PlotPoints  ->  20]; 
h  =  3  10*-7; 

H  =  310*-5; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=300nm  andh=3nm", 

AxesLabel  ->  {"\[Lambda]'*,  'log  \[Sigma]",  ”\[Alpha]  ''}, 

PlotPoints  ->  20]; 
h=  10A-7; 

H=  10A-5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=100nm  andh=lnm", 

AxesLabel  ->  {''\[Lambda]",  "log  \[Sigma]",  M\[Alpha]  "}, 

PlotPoints  ->  20]; 
h  =  3  10A-7; 

H=10A-5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1,  14},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  H=100nm  and  h=3nm", 

AxesLabel  ->  {"\[Lambda]",  "log  \[Sigma]",  "\[Alpha]  "}, 

PlotPoints  ->  20] 
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Program  Darparod2a 
c  =  3  10^10; 

lamb  =  lambmicron/lOM; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf  =  4.2  lO^/rsaO; 

A  =  H/h; 
x  =  kh; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10A-14; 
mfp  =  vftaul ; 

mhopercm  =  mhopercml/(l  +  mfp/h); 
tau  =  mhopercm  taul/mhopercm  1 ; 
sigO  =  1 0^-9  cA2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
m  =  eps^J; 
z  —  m  x; 

qel  =  (2/x)Re[((m  BesselJ[-l,  z]/BesselJ[0,  z])  BesselJ[0,  x]  - 
BesselJ[-l, 

x])/((m  BesselJ[-l,  z]/BesselJ[0,  z])(BesselJ[0,  x]  + 

I  BesselY[0,  x])  -  (BesselJ[-l,  x]  + 

I  BesselY[-l,  x]))  + 

2Sum[((m  BesselJ[j  - 1,  z]/BesselJ[j,  z])BesselJ[j,  x]  - 
BesselJ[j  - 1, 

x])/((m  BesselJ[j  -  1,  z]/BesselJ[j,  z])(BesselJO,  x]  + 
I  Bessel Y[j,  x])  -  (BesselJ[j  - 1,  x]  + 

I  BesselY[j  - 1,  x])),  {j,  1, 19}]]; 
qe2  =  (2/x)Re[(BesselJ[-l,  z]  BesselJ[0,  x]/(m  BesselJ[0,  z])  - 
BesselJ[-l, 

x])/(  (BesselJ[-l,  z]/(m  BesselJ[0,  z]))(BesselJ[0,  x]  + 

I  Bessel  Y[0,  x])  -  (Bessel J[-l ,  x]  + 

I  Bessel Y[-l,  x]»  + 

2Sum[((  BesselJ[j  -  1,  z]/(m  BesselJ[j,  z])  -  j  /(m  z)  + 
j/x)BesselJ[j,  x]  - 
BesselJ[j  - 1, 

x])/(  (BesselJ[j  -  1,  z]/(m  Bessel J[j,  z])  -  j/(m  z)  + 
j/x)(BesselJ[j,  x]  +  I  BesselY[j,  x])  -  (BesselJf 
j  - 1,  x]  +  I  BesselYD  -  1,  x])),  (j,  1, 19}]]; 
qe  =  (qel  +  qe2)/3; 
alpha  1  =  qe/(Pi  (h  +  hO)  10M/4); 
ee  =  (1  -  l/AA2y\5; 

LI  =  (l/eeA2  - 1)(-1  +  Log[(l  +  ee)/(l  -  ee)]/(2ee)); 

L2  =  (l  -  Ll)/2; 

L3  =  L2; 

pi  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kA4(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi)y((Pi/ 
4)H(h+hO)A2  10A4); 
h  =  KTxx; 

mhopercm  1  =  6  lO^  ; 

H  =  4  10A-4; 
h0  =  0; 
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gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=6xl 0*5,H=4\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "Log  h",  n\[Alpha]  PlotPoints  ->  20]; 
mhopercml  =  10*5 ; 
gl  =  Plot3D[ 

Iffalphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*5,H=4\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "Logh",  M\[Alpha]  "},  PlotPoints ->  20]; 
mhopercml  =  10A4 ; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*4,H=4\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "Logh",  "\[Alpha] PlotPoints -> 20]; 
mhopercml  =  10*3 ; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod\[Sigma]=10A3,H=4\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "Logh",  ”\[Alpha]  "},  PlotPoints -> 20]; 
mhopercml  =  6  10A5  ; 

H=10A-4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=6xl0*5,H=l\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "Log  h",  "\[Alpha]  ”},  PlotPoints  ->  20]; 
mhopercml  =  10*5 ; 
gl  =  Plot3D{ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*5,H=l\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "Log  h",  "\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10*4 ; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*4,H=l\[Mu]", 

AxesLabel  ->  {"\[Lambda]",  "Log  h",  "\[  Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10*3 ; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*3,H=l\[Mu]", 

AxesLabel  ->  {"\[Lambda]M,  "Log  h",  "\[  Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  6  10*5; 

H  =  3  10*-5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=6xl0*5JH=300nm", 

AxesLabel  ->  {"\[Lambda]",  "Logh",  "\[  Alpha]  "},  PlotPoints -> 20]; 
mhopercml  =  10*5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*5,H=300nm", 

AxesLabel  ->  {"\[Lambda]",  "Log  h",  "\tAlpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10*4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*4,H=300nm", 

AxesLabel  ->  {"\[Lambda]",  "Logh",  "\[Alpha]  "},  PlotPoints ->  20]; 
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mhopercml  =  10*3; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*3,H=300nm", 

AxesLabel  ->  {^[Lambda]",  "Logh",  "\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  6  10*5; 

H=  1.00001  10*-5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=6xl 0*5,H=  lOOnm", 

AxesLabel  ->  {"\[Lambda]”,  "Logh",  "\[  Alpha]  "},  PlotPoints  -> 20]; 
mhopercml  =  10*5; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*5,H=100nm", 

AxesLabel  ->  {"\[Lambda]M,  "Logh",  "\[ Alpha] "},  PlotPoints -> 20]; 
mhopercml  =  10*4; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*4,H=100nm", 

AxesLabel  ->  {”\[Lambda]",  "Log  h",  "\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10*3; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -5,  -8}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10*3,H=100nm", 

AxesLabel  ->  {"\[Lambda]",  "Log  h",  "\[Alpha]  "},  PlotPoints  ->  20] 


Program  Darparod3a 
c  =  3  10*10; 

lamb  =  lambmicron/1 0*4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10*8/rsa0; 

A  =  H/h; 
x  =  kh; 

taul  =  .22*mhopercml*10*-6*rsa0*3*10*-14; 
mfp  =  vf'taul; 

mhopercm  =  mhopercml/(l  +  mfp/h); 
tau  =  mhopercm  taul/mhopercm  1 ; 
sigO  =  10*-9  c*2  mhopercm; 
sig  =  sigO/  (1  -  I  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
m  =  eps*0.5; 
z  =  mx; 

qel  =  (2/x)Re[((m  BesselJ[-l,  z]/BesselJ[0,  z])  Bessel J[0,  x]  - 
BesselJ[-l, 

x])/((m  BesselJ[-l,  z]/BesselJ[0,  z])(BesselJ[0,  x]  + 

I  BesselYfO,  x])  -  (BesselJ[-l ,  x]  + 

I  Bessel Y[-l,x]))  + 

2Sum[((m  BesselJ[j  -  1,  z]/BesselJ[j,  z])BesselJ[j,  x]  - 
Bessel J[j  - 1, 

x])/((m  Bessel J[j  -  1,  z]/BesselJ[j,  z])(BesselJ[j,  x]  + 
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I  BesselY[j,  x])  -  (BesselJ[j  - 1,  x]  + 

I  BesselYQ  - 1,  x])),  {j,  1, 19}]]; 
qe2  =  (2/x)Re[(BesselJ[-l,  z]  BesselJ[0,  x]/(m  Bessel  J[0,  z])  - 
BesselJ[-l, 

x])/(  (BesselJ[-l,  z]/(tn  Bessel J[0,  z]))(BesselJ[0,  x]  + 

I  BesselY[0,  x])  -  (BesselJ[-l,  x]  + 

I  Bessel Y[-l,x]))  + 

2Sum[((  BesselJ[j  -  1,  z]/(m  BesselJ[j,  z])  -  j  /(m  z)  + 
j/x)BesselJ[j,  x]  - 
Bessel J[j  - 1, 

x])/(  (BesselJ[j  -  1,  z]/(m  Bessel  J[j,  z])  -j/(m  z)  + 
j/x)(BesselJ[j,  x]  + 1  BesselYQ,  x])  -  (BesselJ[ 
j  - 1,  x]  + 1  BesselYQ  - 1,  x])),  {j,  1, 19}]]; 
qe  =  (qel  +  qe2)/3; 
alphal  =  qe/(Pi  (h  +  hO)  10M/4); 
ee  =  (l  -  1/AA2)^.5; 

LI  =  (l/eeA2  -  1X*1  +  Log[(l  +  ee)/(l  -  ee)]/(2ee)); 

L2  =  (l  -Ll)/2; 

L3  =  L2; 

pi  =  (Pi/2)  H hA2(eps  -  l)/(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  H  hA2(eps  -  l)/(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  H  hA2(eps  - 1)/(3  +  3L3(eps  - 1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

k^Absfrl]^  +  Abs[p2]A2  +  Abs[p3]A2)/(l  8  Pi))/((Pi/ 

4)H(h  +  h0)A210A4); 
h  =  10^7; 

mhopercml  =  6  1 0^5 ; 

H=  1.00001  10Axx; 
h0  =  0; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=6xlOA5,h=lmn", 

AxesLabel  ->  {"\[Lambda]w,  "Log  H",  "\[Alpha] "},  PlotPoints  ->  30]; 
mhopercml  =  10^ ; 
gl  =  Plot3D[ 

Iffalphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A5di=lnm”, 

AxesLabel  ->  {"\[Lambda]K,  "Log  H",  ”\[Alpha]  "},  PlotPoints  ->  30]; 
mhopercml  =  10^4 ; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A4di=lnm", 

AxesLabel  ->  {"\[Lambda]",  "Log  H”,  "\[ Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10A3 ; 
gl  =  Plot3D[ 

Iffalphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A3Ji=lnm", 

AxesLabel  ->  {"\[Lambda]",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  6  10*5  ; 
h=  10  10A-7; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=6x10A5,h=10nm", 

AxesLabel  ->  {"\[Lambda]",  "Log  H",  "\[ Alpha] "},  PlotPoints  ->  40]; 
mhopercml  =  10A5 ; 
h=  10  10A-7; 
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gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A5,h=10nm", 

AxesLabel  ->  {"\[Lambda]",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  30]; 
mhopercml  =  10*4 ; 
h  =  10  10A-7; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A43i=l Onm", 

AxesLabel  ->  {^[Lambda]",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10A3 ; 
h  =  10  10^-7; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -6}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A3Jh=10nm", 

AxesLabel  ->  {"\[Lambda]M,  "Log  H",  "\[Alpha] "},  PlotPoints  ->  20]; 
mhopercml  =  6  10^5 ; 
h  =  30  10A-7; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 

14},  {xx,  -3,  -5.5},  PlotLabel  ->  "Metal  Rod  \[Sigma]=6xl0A54i=30nm", 
AxesLabel  ->  {"\[Lambda]M,  "Log  H",  ”\tAlpha]  ”},  PlotPoints  ->  20]; 
mhopercml  =  10A5 ; 
h  =  30  10A-7; 
gl  =Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 

14},  {xx,  -3,  -5.5},  PlotLabel  ->  "Metal  Rod  \[Sigma]=10A5^i=30nm", 
AxesLabel  ->  {"\[Lambda]",  "Log  H",  "\[Alpha]  ”},  PlotPoints  ->  20]; 
mhopercml  =  10*4 ; 
h  =  3010A-7; 
gl  =  Plot3D[ 

If[alphal  <=alpha2,  alphal,  alpha2],  {lambmicron,  1, 

14},  {xx,  -3,  -5.5},  PlotLabel  ->  "Metal  Rod  \[Sigma]=10A4}1=30nm", 
AxesLabel  ->  {"\[Lambda]",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  20]; 
mhopercml  =  10A3  ; 
h  =  3010A-7; 
gl  =  Plot3D[ 

If[alphal  <=  alpha2,  alphal,  alpha2],  {lambmicron,  1, 14},  {xx,  -3,  -5.5}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A33i=30nm", 

AxesLabel  ->  {"\[Lambda]",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  20] 


Program  Darparod4a 
c  =  3  10A10; 

lamb  =  lambmicron/1 0A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10^8/rsa0; 

A  =  H/h; 
x  =  kh; 

taul  =  .22*mhopercml  *10A-6*rsaOA3*10A-14; 
mfp  =  vf*taul; 

mhopercm  =  mhopercml/(l  +  mfp/h); 
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tau  =  mhopercm  taul/mhopercm  1 ; 
sigO  =  lO'M)  c^l  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
m  =  eps'Hl.S; 
z  =  mx; 

qel  =  (2/x)Re[((m  BesselJ[-l,  z]/BesselJ[0,  z])  BesselJfO,  x]  - 
BesselJ[-l, 

x])/((m  BesselJ[-l,  z]/BesselJ[0,  z])(BesselJ[0,  x]  + 

I  BesselYTO,  x])  -  (BesselJ[-l ,  x]  + 

I  Bessel  Y[-l,  x]))  + 

2Sum[((m  BesselJ[j  -  1,  z]/BesselJ[j,  z])BesselJ[j,  x]  - 
Bessel J[j  - 1, 

x])/((m  BesselJU  - 1,  z]/BesselJ[j,  z])(BesselJ[j,  x]  + 

I  BesselY[j,  x])  -  (BesselJ[j  -  1,  x]  + 

I  BesselY[j  -  1,  x])),  (j,  1, 19}]]; 
qe2  =  (2/x)Re[(BesselJ[-l,  z]  BesselJfO,  x]/(m  Bessel J[0,  z])  - 
BesselJ[-l, 

x])/(  (BesselJ[-l,  z]/(m  BesselJ[0,  z])XBesselJ[0,  x]  + 

I  BesselY[0,  x])  -  (BesselJ[-l,  x]  + 

IBesselY[-l,x]))  + 

2Sum[((  BesselJU  -  1,  z]/(m  Bessel  J[j,  z])  -  j  /(m  z)  + 
j/x)BesselJ[j,  x]  - 
BesselJ[j  - 1, 

x])/(  (BesselJU  -  1,  z]/(m  BesselJU,  z])  -  j/(m  z)  + 
j/x)(BesselJU,  x]  + 1  BesselYU,  x])  -  (BesselJ[ 
j  -  1,  x]  +  I  BesselYU  -  1,  x])),  U,  1, 19}]]; 
qe  =  (qel  +  qe2)/3; 
alphal  =  qe/(Pi  (h  +  hO)  10M/4); 
ee  =  (l-l/AA2)A5; 

LI  =  (l/eeA2  -  1X-1  +  Log[(l  +  ee)/(l  -  ee)]/(2ee)); 

L2  =  (1  -  Ll)/2; 

L3  =  L2; 

pi  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Imjpl  +  p2  +  p3]/3  + 

kM(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi)y((Pi/ 

4)H(h  +  h0)A210A4); 
mhopercm  1  =  KKxx ; 
lambmicron  =  14; 
h  =  10Ayy; 
h0  =  0; 

H  =  20  10A-4; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  (yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=14\[Mu]  and  H=20\[Mu]", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]”,  "\[Alpha]"},  PlotPoints  ->  20]; 
H  =  4 10A-4; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5),  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=14\[Mu]  and  H=4\[Mu]", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 
H  =  10A-4; 

gl  =  Plot3Dpf[alphal  <=  alpha2,  alphal ,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=14\[Mu]  and  H=l\[Mu]'\ 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 
H  =  3  10A-5; 
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gl  =  Plot3D[If[alphal  <=  alpha2,  alpha  1,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=14\[Mu]  and  H=300nm", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]"},  PlotPoints  ->  20]; 
H=  1.00001  10^5; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal ,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=l  4\[Mu]  and  H=  1  OOnrn", 
AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[  Alpha]”},  PlotPoints  ->  20] 


Program  Darparod6a 
c  =  3  lCK^lO; 

lamb  =  lambmicron/lOM; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  1 0A8/rsa0; 

A  =  H/h; 
x  =  kh; 

taul  =  .22*mhopercml  *10A-6*rsa0A3*10A-14; 
mfp  =  vftaul; 

mhopercm  =  mhopercml/(l  +  mfp/h); 
tau  =  mhopercm  taul/mhopercml; 
sigO  =  10A-9  0*2  mhopercm; 
sig  =  sigO/  (1  -  I  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
m  =  eps^O.S; 
z  =  mx; 

qel  =  (2/x)Re[((m  Bessel J[-l,  z]/BesselJ[0,  z])  Bessel J[0,  x]  - 
BesselJ[-l, 

x])/((m  BesseU[-l,  z]/BesselJ[0,  z])(BesselJ[0,  x]  + 

I  BesselY[0,  x])  -  (BesselJ[-l,  x]  + 1  BesselY[-l,  x]))  + 
2Sum[((m  Bessel J[j  - 1,  z]/BesselJ[j,  z])BesselJ[j,  x]  - 
Bessel J[j  -  1, 

x])/((m  Bessel J[j  -  1,  z]/BesselJ[j,  z]XBesselJ[j,  x]  + 

I  BesselY[j,  x])  -  (BesselJ[j  -  1,  x]  + 

I  Bessel Y[j  -  1,  x])),  {j,  1, 19}]];  ' 
qe2  =  (2/x)Re[(BesselJ[-l,  z]  BesselJ[0,  x]/(m  BesselJfO,  z])  - 
BesselJ[-l, 

x])/(  (Bessel J[-l,  z]/(m  BesselJ[0,  z])XBesselJ[0,  x]  + 

I  BesselY[0,  x])  -  (Bessel J[-l,  x]  + 1  Bessel Y[-l,  x]))  + 
2Sum[((  BesselJQ  -  1,  z]/(m  BesselJ[j,  z])  -  j  /(m  z)  +  j/x)BesselJ[ 
j,x]- 

BesselJ[j  - 1, 

x])/(  (Bessel J[j  - 1,  z]/(m  Bessel J[j,  z])  -  j/(m  z)  + 
j/x)(BesselJ[j,  x]  + 1  BesselY[j,  x])  -  (BesselJ[ 
j  -  1,  x]  + 1  BesselYjj  -  1,  x])),  {j,  1, 19}]]; 
qe  =  (qel  +  qe2)/3; 
alphal  =  qe/(Pi  (h  +  hOA2/h)  10^4/4); 
ee  =  (1  -  l/AA2y\5; 

LI  =  (l/eeA2  -  iy-1  +  Log[(l  +  ee)/(l  -  ee)]/(2ee)); 

L2  =  (l  -Ll)/2; 

L3  =  L2; 

pi  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3L2(eps  -  1)); 
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p3  =  (Pi/2)  H  hA2(eps  -  l)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kA4(Abs[pl]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi))/((Pi/4)  H  (hA2  + 
hO^)  10^4); 
mhopercml  =  lO''** ; 
lambmicron  =  3; 
h  =  lO^yy; 
h0  =  0; 

H  =  510M; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=3\[Mu],H=5\[Mu]", 

AxesLabel  ->  {"logh",  "log  \[Sigma]",  "\[  Alpha]"},  PlotPoints  ->  20]; 
lambmicron  =  5; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=5\[Mu] ,H=5\[Mu]", 

AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]”},  PlotPoints  ->  20]; 
lambmicron  =  7; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal, alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=7\[Mu],H=5\[Mu]", 

AxesLabel  ->  {"log  h",  ’log  \[Sigma]",  "\[Alpha]”},  PlotPoints  ->  20]; 
lambmicron  =  14; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  3, 6}, 
PlotLabel  ->  "Metal  Rod  \[Lambda]=14\[Mu],H=5\[Mu]", 

AxesLabel  ->  {"log  h",  "log  \[Sigma]",  "\[Alpha]">,  PlotPoints  ->  20]; 


Program  DiskhvH 
0  =  310^0; 

lamb  =  lambmicron/1 0A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf  =  4.2  lO^/rsaO; 

A  =  H/h; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10A-14; 
mfp  =  vf*taul; 
kap  =  h/mfp; 
mhopercm  = 

mhopercml  (1  -  .75(kap  -  kapA3/12)ExpIntegralEi[-kap]  - 

3(1  -  Exp[-kap])/(8kap)  -  (5/8  +  kap/16  -  kapA2/16)Exp[-kap]) ; 
tau  =  mhopercm  taul /mhopercm  1 ; 
sigO  =  10A-9  c*2  mhopercm; 
sig  =  sigO/  (1  - 1  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
a2  =  (2Pi  h  epsA.5/lambXl  -  Sin[th]A2/eps)A.5; 
zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[th]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5)/epsA.5; 
rr  =  I(zl  A2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zl A2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =NIntegrate[ 

Cos[th](l  -  T)10A-4/(h  +  hO),  {th,  10A-4,  Pi/2  -  10A-4}]/(Pi/2); 
ee  =  (1  -  l/A^S; 
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ge  =  (l/ee^2  -  1)^.5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  LI; 

L3  =  l  -LI  -L2; 

pi  =  (Pi/2)  HA2  (h  +  hO)(eps  -  iy(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  HA2  (h  +  hO)(eps  -  iy(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  HA2  (h  +  hO)(eps  -  iy(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

kA4(Abs[p  1  ]A2  +  Abs[p2]A2  +  Abs[p3]A2)/(  1 8  Pi)y((Pi/4)  HA2(h  + 
hO)  10A4); 

(♦lambmicron  =  KK'ss;*) 
h  =  lCK'yy; 

H=  1.00001  lO'bcx; 
h0  =  0; 

mhopercml  =  10^5 ; 
lambmicron  =  3; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=10A5,\[Lambda]=3\[Mu]", 

AxesLabel  ->  {"Log  h",  "Log  H",  "\[Alpha] "},  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  ->  "Metal  Disk  \[Sigma]=l 0A5,\[Lambda]=14\[Mu]", 

AxesLabel  ->  {"Log  h",  "Log  H",  "\[Alpha]  PlotPoints  ->  30]; 
lambmicron  =  3; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

"Metal  Disk(\[Lambda]<=3H  transition)\[Sigma]=10A5,\[Lambda]=3\[Mu]", 
AxesLabel  ->  {"Log  h”,  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

"Metal  Disk(\[Lambda]<=3H  transition)\[Sigma]=10A5,\[Lambda]=14\[Mu]", 
AxesLabel  ->  {"Logh",  "LogH",  "\[Alpha] "},  PlotPoints  ->  30]; 
lambmicron  =  3; 

gl  =  Plot3D[l/(l/alphal  +  l/alpha2),  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

"Metal  Disk(  1  AfAlpha]  1 + 1  A[Alpha]2  \ 
transition)\[Sigma]= 1 0A5,\[Lambda]=3\[Mu] ", 

AxesLabel  ->  {"Logh",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  30]; 
lambmicron  =14; 

gl  =  Plot3D[l/(l/alphal  +  l/alpha2),  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

"Metal  Disk(  1  A[Alpha]  1 + 1  A[Alpha]2  \ 
transition)\[Sigma]=l  0A5,\[Lambda]=  1 4\[Mu] ", 

AxesLabel  ->  {"Logh",  "Log H",  ”\[Alpha] "},  PlotPoints  ->  30]; 


Program  DiskhvH2 
c  =  3  1(^10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf  =  4.2  1 0A8/rsa0; 
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A  =  H/h; 

taul  =  .22*mhopercm  1  *  1 0A-6*rsa0A3  *  1 0*- 1 4; 

mfp  =  vftaul; 

mhopercm  =  mhopercml ; 

tau  =  mhopercm  taul/mhopercm  1 ; 

sigO  =  10^-9  q/^2  mhopercm; 

sig  =  sigO /  (1  - 1  omeg  tau); 

eps  =  1  + 1 4  Pi  sig/omeg; 

a2  =  (2Pi  h  epsA.5/lamb)(l  -  Sin[th]A2/eps)A.  5; 

zl  =  (l/Cos[th]  +  Cos[th])/2; 

z2  =  (1/(1  -  Sin[fli]A2/eps)A.5  +  (1  -  Sin[th]A2/eps)A.5yepsA  5; 
rr  =  I(zl  A2  -  z2A2)Tan[a2]/(2zl  z2  -  I(zlA2  +  z2A2)Tan[a2]); 
tt  =  2zl  z2/(2zl  z2  Cos[a2]  -  I(zl A2  +  z2A2)Sin[a2]); 

R  =  Abs[rr]A2; 

T  =  Abs[tt]A2; 
alphal  =  NIntegrate[ 

Cos[th](l  -  T)10A-4/(h  +  hO),  {th,  10M,  Pi/2  -  10A-4}]/(Pi/2); 
ee  =  (l-l/AA2)A.5; 
ge  =  (l/eeA2  - iy\ 5; 

LI  =  ge/(2eeA2)(Pi/2  -  ArcTan[ge])  -  geA2/2; 

L2  =  LI; 

L3  =  1  -  LI  -  L2; 

pi  =  (Pi/2)  HA2  (h  +  hOXeps  -  iy(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  HA2  (h  +  hOXeps  -  iy(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  HA2  (h  +  hOXeps  -  iy(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

k^Abstpl]^  +  Abs[p2]A2  +  Abs(p3]A2y(18  Pi)y((Pi/4)  HA2(h  + 
hO)  10A4); 

(♦lambmicron  =  lO^s;*) 
h=  10Ayy; 

H=  1.00001  lO'Tcx; 
h0  =  0; 

mhopercml  =  10*5 ; 
lambmicron  =  3; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  -> 

"Metal  Disk(No  size  effects)\[Sigma]=10A5,\[Lambda]=3\[Mu]", 
AxesLabel  ->  {"Log  h",  "Log  H",  "\[Alpha]  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal, alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  -> 

"Metal  Disk(No  size  effects)\[Sigma]=10A5,\[Lamb<la]=14\[Mu]", 
AxesLabel  ->  {"Logh",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  30]; 
lambmicron  =  3; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  -> 

"Metal  Disk(No  size  effects)\[Sigma]=10A5,\[Lambda]=3\[Mu]", 
AxesLabel  ->  {"Logh",  "Log  H",  ”\[Alpha]  "},  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  -> 

"Metal  Disk(No  size  efFects)\[Sigma]=10A5,\[Lambda]=14\[Mu]", 
AxesLabel  ->  {"Logh",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  30]; 
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Program  RodhvH2 
c  =  3  10^10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf=4.2  10*8/rsa0; 

A  =  H/h; 
x  =  kh; 

taul  -  .22*mhopercml*10*-6*rsa0*3*10*-14; 
mfp  =  vPtaul; 

mhopercm  =  mhoperctnl/(l  +  mfp/h); 
tau  =  mhopercm  taul /mhopercm  1; 
sigO  =  10**9  c*2  mhopercm; 
sig  =  sigO/  (1  -  I  omeg  tau); 
eps  =  1  + 1 4  Pi  sig/omeg; 
m  =  eps*0.5; 
z  =  m  x; 

qel  =  (2/x)Re[((m  BesselJ[-l ,  z]/BesselJ[0,  z])  BesselJ[0,  x]  - 
Bessel  J[-l, 

x])/((m  Bessel J[-l,  z]/BesselJ[0,  z])(BesselJ[0,  x]  + 

I  Bessel Y[0,  x])  -  (BesselJ[-l,  x]  + 

I  Bessel  Y[-l,  x]))  + 

2Sum[((m  BesselJ[j  - 1,  z]/BesselJ[j,  z])BesselJO,  x]  - 
Bessel  JO  - 1, 

x])/((m  BesselJ[j  -  1,  z]/BesselJ[j,  z])(BesselJO,  x]  + 
I  Bessel  Y[j,  x])  -  (BesselJO  -  1,  x]  + 

I  BesselYO -l,x])),{j,  1,19}]]; 
qe2  =  (2/x)Re[(BesselJ[-l,  z]  Bessel J[0,  x]/(m  Bessel J[0,  z])  - 
Bessel  J[-l, 

x])/(  (BesselJ[-l,  z]/(m  BesselJ[0,  z]))(BesselJ[0,  x]  + 

I  Bessel Y[0,  x])  -  (BesselJ[-l,  x]  + 

I  BesselY[-l,  x]»  + 

2Sum[((  BesselJ[j  - 1,  z]/(m  BesselJO,  z])  -  j  /(m  z)  + 
j/x)BesselJ[j,  x]  - 
BesselJO  -  1, 

x])/(  (BesselJ[j  -  1,  z]/(m  BesselJO,  z])  -  j/(m  z)  + 
j/x)(BesselJO,  x]  +  I  BesselYO,  x])  -  (BesselJ[ 
j  -  1,  x]  + 1  BesselYO  - 1,  x]»,  (j,  1, 19}]]; 
(•combine  results  for  random  orientation  in  plane  with  G  = 

S/4  =  2 A/4  =  A/2  for  3D  random  orientation*) 
qe  =  (qel  +  qe2)/4; 
alphal  =  qe/(Pi  (h  +  hO)  10*4/4); 
ee  =  (1  -  1/A*2)*.5; 

LI  =  (l/ee*2  -  1X-1  +  Log[(l  +  ee)/(l  -  ee)]/(2ee)); 

L2  =  (l  -  Ll)/2; 

L3  =  L2; 

pi  =  (Pi/2)  H  h*2(eps  -  iy(3  +  3Ll(eps  - 1)); 
p2  =  (Pi/2)  H  h*2(eps  -  l)/(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  H  h*2(eps  - 1)/(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Imfpl  +  p2  +  p3]/3  + 

k*4(Abs[p  1  ]*2  +  Abs[p2]*2  +  Abs[p3]*2)/(1 8  Pi))/((Pi/ 

4)  H  (h  +  h0)*2  10*4); 
h  = 10*yy; 

H=  1.00001  10*xx; 
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h0  =  0; 

mhopercml  =  10^5 ; 
lambmicron  =  3; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5),  {xx,  -1,  -5}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=  1  O^AfLambdaHVIMu] ", 

AxesLabel  ->  {"Logh",  "Log  H",  "\[Alpha]  PlotPoints  ->  30]; 
lambmicron  =14; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  ->  "Metal  Rod  \[Sigma]=10A5,\[Lambda]=14\[Mu]", 

AxesLabel  ->  {"Log  h",  "Log  H”,  "\[Alpha]  ">,  PlotPoints  ->  30]; 
lambmicron  =  3; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

"Metal  Rod(\[Lambda]<=3H  transition)\tSigma]=10A5,\[Lambda]=3\[Mu]", 
AxesLabel  ->  {"Log  h",  "Log  H",  "\[Alpha]  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal ,  alpha2],  {yy,  -8,  -5},  {xx,  -1,-5}, 

PlotLabel  -> 

"Metal  Rod(\[Lambda]<=3H  transition)\tSigma]=10A5,\[Lambda]=14\[Mu]", 
AxesLabel  ->  {"Log  h",  "Log  H",  "\{Alpha]  '*},  PlotPoints  ->  30]; 
lambmicron  =  3; 

gl  =  Plot3D[l/(l/alphal  +  l/alpha2),  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

’Metal  Rod(lA[Alpha]l+lA[Alpha]2  \ 
transition)\[Sigma]=10A5,\{Lambda]=3\[Mu]", 

AxesLabel  ->  {"Logh",  "Log  H",  "\[Alpha] PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =Plot3D[l/(l/alphal  +  l/alpha2),  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

"Metal  Rod(l A[ Alpha]  1+1  A[Alpha]2  \ 
transition)\[Sigma]=l  0A5,\[Lambda]=  1 4\[Mu] ", 

AxesLabel  ->  {"Log  h”,  "Log  H",  ”\[Alpha] PlotPoints  ->  30]; 


Program  RodhvH3 
c  =  3  10*10; 

lamb  =  lambmicron/10A4; 
omeg  =  2  Pi  c/lamb; 
k  =  2  Pi  /  lamb; 
rsaO  =  2.5; 
vf  =  4.2  10A8/rsaO; 

A  =  H/h; 
x  =  kh; 

taul  =  .22*mhopercml*10A-6*rsa0A3*10A-14; 

mfp  =  vftaul ; 

mhopercm  =  mhopercml; 

tau  =  mhopercm  taul /mhopercm  1 ; 

sigO  =  10A-9  c*2  mhopercm; 

sig  =  sigO/  (1  -  I  omeg  tau); 

eps  =  1  + 1 4  Pi  sig/omeg; 

m  =  eps*0.5; 

z  =  mx; 

qel  =  (2/x)Re[((m  BesselJ[-l,  z]/BesselJ[0,  z])  BesselJfO,  x]  - 
Bessel  J[-l, 
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x])/((m  BesselJ[-l,  z]/BesselJ[0,  z])(BesselJ[0,  x]  + 

I  BesselY[0,  x])  -  (BesseIJ[-l,  x]  + 

I  Bessel Y[-l,  x]))  + 

2Sum[((m  Bessel  JO  - 1,  z]/BesselJO,  z])BesselJO,  x]  - 
BesselJ[j  -  1, 

x])/((m  Bessel  JO  -  1,  z]/BesselJO,  z])(BesselJO,  x]  + 

I  Bessel  YO,  x])  -  (Bessel JO  - 1,  x]  + 
IBesselY0-l,x])),G,l,19}]]; 
qe2  =  (2/x)Re[(BesselJ[-l,  z]  Bessel J[0,  x]/(m  Bessel J[0,  z])  - 
BesselJ[-l, 

x])/(  (Bessel  J[-l,  z]/(m  Bessel  J[0,  z])XBesselJ[0,  x]  + 

I  Bessel  Y[0,  x])  -  (BesselJ[-l,  x]  + 

I  BesselY[-l,  x]))  + 

2Sum[((  Bessel  JO  -  1,  z]/(m  Bessel  JO,  z])  -  j  /(m  z)  + 
j/x)BesselJO,  x]  - 
Bessel  J[j  - 1, 

x])/(  (Bessel JO  - 1,  z]/(m  BesselJO,  z])  -  j/(m  z)  + 
j/xXBesselJO,  x]  + 1  BesselYO,  x])  -  (BesselJ[ 
j  -  1,  x]  + 1  BesselYQ  -  1,  x])),  0, 1, 19}]]; 

(^combine  results  for  random  orientation  in  plane  with  G  = 

S/4  =  2A/4  =  A/2  for  3D  random  orientation*) 
qe  =  (qel  +  qe2)/4; 
alphal  =  qe/(Pi  (h  +  hO)  10M/4); 
ee  =  (l-l/AA2)A.5; 

LI  =  (l/eeA2  -  iy-1  +  Log[(l  +  ee)/(l  -  ee)]/(2ee)); 

L2  =  (l  -Ll)/2; 

L3  =  L2; 

pi  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3Ll(eps  -  1)); 
p2  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3L2(eps  -  1)); 
p3  =  (Pi/2)  H  hA2(eps  -  iy(3  +  3L3(eps  -  1)); 
alpha2  =  (k  Im[pl  +  p2  +  p3]/3  + 

k^Absfpl]^  +  Abs[p2]A2  +  Abs[p3]A2)/(18  Pi)y((Pi/ 

4)H(h  +  hO)A2  10A4); 
h  =  10Ayy; 

H=  1.00001  10Axx; 
h0  =  0; 

mhopercml  =  10A5 ; 
lambmicron  =  3; 

gl  =Plot3D[If[alphal  <=  alpha2,  alphal ,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 
PlotLabel  ->  "Metal  Rod(No  size  effects)\[Sigma]=10A5,\[Lambda]=3\[Mu]", 
AxesLabel  ->  ("Log  h",  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[If[alphal  <=  alpha2,  alphal,  alpha2],  (yy,  -8,  -5),  {xx,  -1,  -5}, 
PlotLabel  -> 

"Metal  Rod(No  size  effects)\[Sigma]=10A5,\[Lambda]=14\[Mu]", 

AxesLabel  ->  ("Logh",  "Log  H",  "\[Alpha]  PlotPoints  ->  30]; 
lambmicron  =  3; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  ->  "Metal  Rod(No  size  effects)\[Sigma]=10A5,\[Lambda]=3\[Mu]", 
AxesLabel  ->  {"Log  h”,  "Log  H",  "\[Alpha]  "},  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[If[lamb  <=  3H,  alphal,  alpha2],  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  -> 

"Metal  Rod(No  size  efFects)\[Sigma]=10A5,\[Lambda]=14\[Mu]", 

AxesLabel  ->  {"Log  h",  "Log  H",  M\[Alpha]  "},  PlotPoints  ->  30]; 
lambmicron  =  3; 
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gl  =  Plot3D[l/(l/alphal  +  l/alpha2),  {yy,  -8,  -5},  {xx,  -1,  -5}, 

PlotLabel  ->  "Metal  Rod(No  size  effects)\[Sigma]=10A5,\[Lambda]=3\[Mu]", 
AxesLabel  ->  {"Log  h",  "Log  H",  "\[Alpha] "},  PlotPoints  ->  30]; 
lambmicron  =  14; 

gl  =  Plot3D[l/(l/alphal  +  l/alpha2),  {yy,  -8,  -5},  {xx,  -1,  -5), 

PlotLabel  -> 

"Metal  Rod(No  size  effects)\[Sigma]=10A'5,\[Lambda]=14\[Mu]", 

AxesLabel  ->  {"Log  h”,  "Log  H”,  "\[Alpha] "},  PlotPoints  ->  30]; 
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